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1 


The detections of gravitational waves (GW) by LIGO/Virgo collaborations provide 
various possibilities to physics and astronomy. We are quite sure that GW observations 
will develop a lot both in precision and in number owing to the continuous works for 
the improvement of detectors, including the expectation to the newly joined detector, 
KAGRA, and the planned detector, LIGO-India. In this occasion, we review the funda- 
mental outcomes and prospects of gravitational wave physics and astronomy. We survey 
the development focusing on representative sources of gravitational waves: binary black 
holes, binary neutron stars, and supernovae. We also summarize the role of gravitational 
wave observations as a probe of new physics. 
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1. Introduction 


Table 1 Observing terms of LIGO/Virgo detectors. Detectors column: H1 at Hanford, 
Washington, and L1 at Livingston, Louisiana. V1 at Cascina, Italy. 


Observation Period detectors events reported 
O1 Sep 2015 to Jan 2016 | HI, L1 3 BBHs 
O2 Nov 2016 to Sep 2017 | H1, L1 3 BBHs 
O2 Aug 2017 H1, L1, V1 4 BBHs, 1 BNS 
O3a Apr 2019 to Sep 2019 | H1, L1, V1 many BBHs, 1 BNS 
O3b Nov 2019 to Apr 2020 | H1, L1, V1 N/A 


1.1. Gravitational wave from binary black holes 


The first gravitational wave (GW) event from the binary black holes (BBHs), GW150914, 
which was discovered by LIGO/Virgo collaboration brought us many surprises [1]. Regarding 
binary neutron star (BNS) coalescences, it was believed that there is not much uncertainty in 
estimating their event rate compared with the case of BBHs. Even for BNS coalescences, the 
uncertainty in the estimated event rate was expected to be about one order of magnitude. 
On the other hand, although it is conceivable that BBHs should exist in our universe, the 
estimation of the GW event rate had no observational support, differently from the case of 
BNSs. Hence, BBHs were not considered to be a promised GW source. It was somewhat 
surprising that such a system was the first directly detected GW source. However, it was 
even more surprising that those black holes (BHs) were heavier than many people expected, 
i.e., each component BH has a mass about 30M. The mass of the BH candidates found so 
far by electromagnetic observations is estimated to be at most about 20M; [2]. Although 
there were some candidate objects that did not contradict with even heavier mass within the 
observational error, there is no object that can be said to have about 30M. The observation 
of GWs has changed our understanding about the mass distribution of BHs in the universe |3, 
4]. Since then, the number of BBH observations officially announced by the LIGO/Virgo 
collaborations has been increasing [5]. Not all of them are massive BHs, but the succeeding 
observations indicate that the first GW event, GW150914, was not so special. 

As a result of this discovery, how those BBHs detected by GWs were formed in the history 
of the universe has become one of the great mysteries of astrophysics. A variety of BBH 
formation scenarios have been proposed [6]. The most promising candidate for the forma- 
tion channel would be the one through the binary interaction between heavy star binaries 
formed from gases with a low metallicity. Even for this leading scenario, there are various 
uncertainties: the binary formation rate from the gas cloud, stellar evolution including mass 
loss processes, binary interaction to form close binaries, formation process of BHs, and so on. 
The binary evolution can be so complex, changing its orbital radii and eccentricity. There- 
fore, it is a theoretical challenge to predict the distributions of BH masses, orbital separation 
and eccentricity at the formation of a BH binary, to make comparison with observations. 
As alternative scenarios, we can think of the formation of BH binaries in star clusters [7] 
and their formation in the primordial universe [8]. Many variants of these scenarios can be 
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thought of. In the future, as the details of the theoretical model become better understood, 
and as GW observations accumulate the information about the event rate, the mass and 
distance distributions, the distribution of BH spins, and so on, we would be able to identify 
the contributions from respective formation channels of BBHs. 


1.2. Gravitational wave from binary neutron stars 


The first GW event from a BNS merger GW170817 observed on August 17, 2017, marked 
a splendid opening of GW astronomy [9]. BNS coalescence was counted as a guaranteed 
GW source. The discovery was followed by simultaneous observation of gamma-ray bursts, 
identification of the electromagnetic counterpart by optical and infrared telescopes [10], long- 
term follow-up observations by radio waves and X-rays, detection of apparent superluminal 
motion [11, 12]. How BNS coalescence proceed was revealed to be something close to what 
was theoretically predicted. In the study of BNS coalescence, weak gravity approximation 
is not appropriate any more. Hence, numerical relativity is an indispensable research tool 
with the knowledge of equation of state for the high density nuclear matter. At the same 
time various radiation fields, photons and neutrinos, originating from high density matter 
must be solved and the magnetic fields also can play an important role in the simulations. 
Therefore combining various techniques for numerical simulation is required. In the case of 
GW170817, from its optical and infrared light curves, it is thought that the ejecta mass was 
about 0.05 M and a part of them possessed high neutron excess, resulting in production of 
significant amount of r-process elements. The amount of the produced r-process elements 
may be sufficient to explain almost all of them in the universe. In addition, numerical sim- 
ulations show that the merged binary becomes a massive neutron star (NS) supported by 
rotation, and the material stripped around it forms a disk. The simultaneous (1.7 second 
after the coalescence time) observation of a short gamma-ray burst (SGRB) indicates the 
connection between SGRB and BNS coalescence. Although the observed SGRB was extraor- 
dinary weak, this can be understood because the SGRB jet was not pointing at the line of 
sight toward us [13, 14]. This interpretation gained a strong support because apparent super- 
luminal motion was observed by radio waves [11, 12]. The detection of apparent superluminal 
motion is also a strong evidence for that a relativistic jet has been launched. In addition, the 
GW observations at this event also provided restrictions on the equation of state of dense 
nuclear matter, mostly from the constraint on the tidal deformability of NSs just before the 
coalescence [9]. Here again, numerical relativity technique played an important role. It is 
also an important theoretical challenge to derive the equation of state that can be compared 
with phenomena from the first principle of QCD. High-energy gamma-rays, neutrinos, and 
cosmic rays were not detected in GW170817 [15-18], although BNS mergers and SGRBs 
are considered to be sources of these particles [19-22]. We do not further discuss these 
multi-messenger signals in this review. 


1.9. Gravitational wave from Supernovae 

Although GWs from supernovae have not been detected yet, they are important targets of 
gravitational wave observations. After more than a half-century of continuing efforts ever 
since the seminal work by Colgate and White (1966) [23], theory and neutrino radiation- 
hydrodynamics simulations of core-collapse supernovae (CCSNe) are now converging to a 
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point that multi-dimensional (multi-D) hydrodynamics instabilities play a key role in the 
neutrino mechanism [24], the most favoured scenario to trigger explosions. In fact, self- 
consistent simulations in three spatial dimensions (3D) now report shock revival of the 
stalled bounce shock into explosion by the multi-D neutrino mechanism (see [25-28] for 
reviews). 

From observational point of view, mounting evidence from multi-wavelength 
electromagnetic-wave (EM) observations (such as blast/ejecta morphologies, line profiles, 
polarizations) are pointing toward CCSNe being generally aspherical (e.g., [29, 30] and ref- 
erences therein). Nevertheless, these EM signals are secondary as a probe into the supernova 
inner-workings, providing information only after a shock break-out from the massive stars. 
On the other hand, the GWs from CCSNe are direct observables, which imprint a live infor- 
mation of the central engine. Furthermore coincident detection of CCSN neutrinos with GWs 
is important not only for significantly enhancing the detectability of the GWs (e.g., [31]), 
but also for unraveling the hydrodynamics evolution from the onset of core-collapse toward 
the forming compact objects (neutron stars and black holes, see [32] for a review). Current 


estimates of CCSN rates in our Milky way predict one event every ~ 40 + 10 year [33]. While 
rare, they provide us unique opportunity to study CCSNe using a full set of multi-messenger 
observables including GWs, neutrinos, and nuclear gamma-rays [34]. 


1.4. Test of general relativity 


The direct detection of GWs has opened a new window for tests of general relativity (GR). 
In order to detect GW signals, we need theoretical prediction for the GW waveform in 
advance. In more complex theory of gravitation that extends GR, theoretically predicting the 
waveform is more difficult. Our current knowledge about the waveform in modified gravity 
is limited, but some possible deviations in the waveform from GR have been predicted. GW 
observations so far have not suggested any deviation from GR [35]. However, the lack of 
deviation itself can have important physical meaning. Besides constraints on the deviation 
in the waveform, an important constraint was also derived from the speed of propagation of 
GWs, i.e., from the detection of short gamma-ray burst associated with GW170817 event 
at a distance of about 40 Mpc. The fractional deviation of the propagation speed of GWs 
from that of light is less than œ~ 107? [36], which excludes various possible extensions of 
GR [87, 38, 40-43]. 

GWs from binaries can be used as the standard siren [44], providing an alternative method 
to measure the cosmic expansion history. GW170817 already gives us an estimate of the 
Hubble constant [45, 46], but we do not further discuss the aspects of testing cosmological 
models using GWs in this review. 


1.5. Organization of this article 

In this review paper we overview the vast area of newly developing gravitational wave physics 
and astronomy mentioned above. This paper is organized as follows in Sec. 2, 3 and 4, we 
focus on three major sources of gravitational waves, i.e., BBHs, BNSs, and SNe. We discuss 
various topics emanating from the detection of gravitational waves, including the future 
prospects. In Sec. 5 we focus more on the new physics examined by using GWs as a probe. 
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We discuss the test of gravity and implication to the dark matter physics. Sec. 6 is devoted 
to a short conclusion. 
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2. Binary Black Holes 

2.1. Gravitational waves from binary black holes 

LIGO/Virgo has reported GWs from BBH mergers as a result of O1 and O2 observing run [5] 
that is called “Gravitational-Wave Transient Catalog" (GWTC)-1. The best fit parameters 
for these binaries are listed in Table 2. We have large errors in the parameter estimation, 
especially for the individual masses and the effective spin parameter, because the signal-to- 
noise ratios (SNRs) are not so high. There are results reported by alternative searches. One 
of them reports other candidates of GW events [47]. 


Table 2 BBH merger events observed during LIGO/Virgo O1 and O2. We show only the 
median values for the individual masses mı and m», chirp mass M, effective spin xeg, final 
BH mass Mg, and its angular momentum parameter ag. These numbers are of source frame. 
The estimated distance in Mpc and red-shift z are also shown. The 90% credible intervals 
have been reported in Table III of Ref. [5]. 


Event mı/Mo ma3/Mo M/Mo Xer | Ma/Mo ag | distance (Mpc,z) 


GW150914 35.6 30.6 28.6 -0.01 63.1 0.69 | 430 Mpc 0.09 
GW151012 23.2 13.6 15.2 0.05 35.7 0.67 | 1060 Mpc 0.20 
GW151226 13.7 7.7 8.9 0.18 20.5 0.74 | 440 Mpc 0.09 
GW170104 30.8 20.0 21.4 -0.04 49.1 0.66 | 960 Mpc 0.19 
GW170608 11.0 7.6 7.9 0.03 17.8 0.69 | 320 Mpc 0.07 


GW170729 50.2 34.0 35.4 0.37 80.3 0.81 | 2750 Mpc 0.48 
GW170809 39.0 23.8 24.9 0.08 56.4 0.70 | 990 Mpc 0.20 


GW170814 30.6 25.2 24.1 0.07 53.4 0.72 | 580 Mpc 0.12 
GW170818 35.4 26.7 26.5 -0.09 59.8 0.67 | 1020 Mpc 0.20 
GW170823 39.5 29.0 29.2 0.09 65.6 0.71 | 1850 Mpe 0.34 


The obtained mass and spin distributions are summarized in Fig. 1. The chirp mass is 
defined by 


M = PM (1) 


with M = m; + m3 and u = mym2/M, where m, and mə are the component masses. The 
effective spin parameter is defined by 


1 S S ^ 
vot = G7 (S42) a. (2) 


M My mg 


where S, and S» are spin angular momentum vectors of the respective bodies and Ê is the 
unit vector pointing at the direction of the orbital angular momentum of the binary. 

The reported event rate by the LIGO/Virgo Collaborations is ranging from 9.7 to 101 
Gpc ?yr-! for BBH mergers. The masses of the observed BHs are relatively large compared 
with the masses of stellar mass BH candidates found by X-ray observations (see, e.g., Ref. [2] 
and references therein). The origin of BBHs is an important new issue of debate raised by 
the GW observations. A promising formation path is the one through the stellar evolution 
of low metallicity stars, which is the main topic of Section 2.2. There are various formation 
channels and the main focus is whether or not the stellar dynamics in the star cluster plays 
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Fig. 1 Distribution of the median values of the chirp mass M, mass ratio q = m2/mj, 
effective spin Xer, and individual masses (which are redundant though) of BBHs observed 
during LIGO/Virgo O1 and O2. It should be noted that there are large errors in the current 
observations. Even for the chirp mass which is the most sensitive parameter in the GW 
waveform of BBH mergers, the 90% credible interval gives < 20% difference in the estimation. 


a crucial role. An alternative is the formation of BH binaries in the early universe, which 
is the so-called primordial BH scenario. The formation process of primordial BH binaries is 
discussed in Section 5.2. 


2.2. Theoretical study on binary black hole formation 


Mergers of BBHs are the dominant population of gravitational wave sources. Currently, two 
evolutionary pathways are envisaged as origins of BBHs. First is the so-called in-situ BBH 
formation scenario, where massive star binaries are formed and each member star becomes a 
black hole (BH) after stellar evolution. Second is the so-called ex-situ scenario, where single 
stars formed in a dense star cluster become BHs after stellar evolution and those BHs become 
binaries as a result of dynamical interaction. Here, we describe our current understanding 
of those scenarios in this order. 


2.8.1. Formation of isolated stellar binaries. 
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Fig. 2 An example of the latest 3D radiation-hydrodynamics simulations of the primordial 
star formation, in which the metallicity is exactly zero (from Sugimura et al. submitted [48]), 
showing a snapshot at the epoch when massive protostars in a binary system accrete the 
gas under the radiative feedback. Presented are the volume-rendered 3D distribution of the 
mass density, together with ionization fronts (yellow surfaces) and protostars (black dots). 


A. Binary star formation at low metallicities 

Whereas massive stellar binaries are possible astrophysical progenitors of BH-BH binaries 
that finally merge, their formation process is still uncertain. In particular, major focus is put 
on formation of such systems at low metallicities Z S 0.1 Zo, which are required for yielding 
massive (~ a few x10 Mo) BHs avoiding significant mass loss during the stellar evolution. 
Observations suggest that a large fraction of massive stars are born in binaries in the solar 
neighborhood. Some of them are in very tight systems with separations of only < 0.1 AU, 
and are possible candidates that finally evolve into BH binaries that merge within the Hubble 
time. On the other hand, since the delay time until the BBH merger is potentially as long as 
the Hubble timescale, depending on the initial separation, a part of the progenitor stars may 
have been born in the early universe. Therefore, our goal is to elucidate the formation process 
of massive and tight binary stars in a variety of low-metallicity environments, including the 
local and early universe. 

A straightforward method which ultimately leads us to the goal is directly following the 
evolution with numerical simulations. One of the key processes is the gas gravitational frag- 
mentation that potentially leads to bearing multiple stars, some of which are to be in binaries. 
Magnetic fields, if any, should be important because it transfers the angular momentum of the 
gas, affecting the orbital evolution of newly forming binaries. The stellar radiative feedback 
is also critical to determine how massive binaries finally appear, halting the mass accretion 
onto the stars [49-51]. We consequently need to solve the interplay between gas dynamics, 
magnetic fields, and radiative fields in full 3D spatially resolving a very wide dynamic range. 
Ultimately integrating the evolution over ~ Myr is required to cover the whole evolution of 
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the star formation process. Although this is still challenging, recent studies are advancing to 
incorporate the above processes as much as possible with limited computational resources. 
For instance, Fig. 2 shows a snapshot in our recent 3D radiation hydrodynamic simulations 
using the adaptive-mesh-refinement code SFUMATO [48]. In this particular case, we follow 
the long-term (~ 0.1 Myr) evolution in the protostellar accretion stage of the primordial 
star formation. We see a multiple stellar system including a massive binary system that has 
appeared after the disk fragmentation in the earlier stage. The double bipolar HII regions 
are growing from the two massive stars that are still accreting the surrounding gas under 
the radiative feedback. Although this provides a glimpse of the binary formation in the early 
universe, the newborn binary has a large separation: the two massive stars exciting the bipo- 
lar HII regions are separated by more than 10? AU. The formation process of the massive 
and tight binaries, some of which should evolve to the gravitational wave sources, is yet to 
be revealed. Obviously further studies and simulations are necessary for explaining how such 
systems are formed. We describe our recent studies and future prospects in what follows. 


B. Role of magnetic fields in close binary formation 

Observations of nearby star-forming regions indicate that a large fraction of stars are 
members of binary systems [52] and the binary frequency is even higher for high-mass stars 
than for low-mass ones [53]. Theoretically, numerical simulations have shown that a binary 
stellar system forms as a result of fragmentation of a molecular cloud core during the grav- 
itational collapse [54]. In rapidly rotating and gravitationally unstable cores, fragmentation 
occurs more easily [55, 56] and the accretion rate is higher after the protostars are formed. 
The rotation velocity increases and the cloud becomes disk-like during the collapse owing to 
the angular momentum conservation. This promotes fragmentation and subsequent binary 
formation. 

According to numerical simulations, the binary separation continues to increase with time 
as a parcel of gas with larger angular momentum falls onto the central region or proto-binary 
system at later times [54, 57] and the outcome is a wide binary system with a separation 
typically greater than >> 100 AU. In contrast, observationally high-mass stars tend to be in 
close-binary systems with the separation 1-10 AU [52]. This discrepancy between theory and 
observations can be circumvented if we take non-ideal magnetohydrodynamics (MHD) effect 
into account. 

Recent non-ideal MHD simulations for present-day star formation showed that the excess 
angular momentum is transported from the proto-binary system by magnetic braking and 
by outflow launching and, as a consequence, the close binary system can be formed in a 
highly unstable magnetized cloud. Figures 3 is a snapshot of a proto-binary system in an 
example of such simulations [58]. In the figure, we can see that powerful jet is driven by each 
of the protostars. The separation of the binary system is about 10 AU at this epoch. On the 
equatorial plane, both the circumstellar and circumbinary disks can be confirmed (the left 
panel of Fig. 3). Two cavity-like structures created by the jets are seen in the middle panel 
of Fig. 3. The velocity of the jets driven by the protostars exceeds > 100kms-!. Those 
high-velocity components are enclosed by the low-velocity outflow of ~ I0 kms"! driven 
by the circumbinary disk. The angular momentum of the binary system is extracted from 
the system both by the large-scale wide-angle outflows and small-scale collimated jets. The 
binary jet/outflow structure is also reported from recent ALMA observations [59]. 
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Fig. 3  Proto-binary formation in the MHD simulation. Shown are density (color) and 
velocity (arrows; middle and right panels) distributions on the z = 0 (left) and y = 0 (middle 
and right) plane at three different spatial scales. The elapsed time after the onset of prestellar 
cloud ¢ and that after protostar formation tp, are described in the left panel. 


C. Stellar system formed from metal-free gas 

Next we see the nature of first stellar systems formed from the primordial gas and discuss 
possible formation of first star binaries. The collapse of a slightly gravitationally unstable 
gas cloud (ng = 1.4 x 10*cm^? and T = 200K) is simulated by the SPH method. Armed 
with elaborate sink particle creation scheme and usage of the barotropic equation of state 
pre-calculated by a one-zone model, we succeeded in extending the time integration up to 
25000 yr in the accretion phase, ~ 17 times longer than the previous calculation [60]. 

Figure 4 shows the number of fragments, i.e., sinks, as functions of the scaled time nor- 
malized by the free-fall time at the density ni, = 10!°cm~?, above which the temperature 
rises adiabatically. The results can be scaled using the scaled time according to the scale-free 
nature of gravity [60]. Note that the integration time of the present simulation is also the 
longest so far in terms of the scaled time. We can see a simple trend that the number of frag- 
ments is proportional to the elapsed time after the central density exceeds ny. This trend 
has already been known previously [60], but now we confirm that the relation holds even 
after 10 times longer integration time, although the simulation is done for one particular 
realization. This result is also generally consistent with results of other simulations in the 
literature. 

Several systems of binaries are identified in this cluster. Figure 5 shows the time evolution 
of their separations. Different colors correspond to different pair of binaries. The members 
of the binaries easily change as a result of the close encounter/merger with other stars. The 
thick red curve at the top of the panel shows the typical size of this cluster, which grows 
in proportion to the elapsed time after the first protostar formation, as is predicted by the 
conservation of angular momentum and the gas distribution at the end of the collapse phase. 
On the other hand, the binary separations also grow when the binary systems acquire the 
angular momentum as a result of the gas accretion, which is, however, counterbalanced by the 
effect of encounters with other stars that makes the binary tighter. Hence, the separations are 
kept rather constant during the computation. The separations observed in this particular 
run are ~ a few 10 AU, which is too wide to end up with BBHs that coalesce within a 
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Fig. 4 Number of fragments in a first-star-forming cloud as function of the scaled time. 
We assume paq is 2 x 107?(g/cm-?), above which the gas behaves adiabatic. Red line with 
small triangles denotes the present result, and the thin lines correspond to our previous 
result |60]. Results in the literature[61-71]| are also shown. 
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Fig. 5 The separations of binaries identified in the present simulation as functions of 
time. Different colors denote different binaries. T'hick red solid curve shows the mean square 
average of the distance between the stars and the center of mass of the stellar system whereas 
the green thin curve shows the distance in proportion to the elapsed time since the first sink 
formation. 
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Hubble time. In the present simulation, however, we adopt a lower threshold density than 
the true adiabatic density for computational reasons and the actual binary separation should 
be scaled from our nominal value. After the scaling, the time should be 100 times shorter 
and the separation ~ 70 times less. Thus, our binary separation corresponds to a few times 
0.1 AU after scaling, although the integrated time is only 250 yr. If this separation is kept 
constant for another ~ 5000 yr when the surrounding gas is swept away by the radiative 
feedback, these binaries could be progenitors of merging BBHs observed in gravitational 
waves. 

To draw firm conclusions about the nature of formed stellar systems, we need to extend the 
numerical simulation in time until the radiative feedback from the protostar terminates the 
accretion to the system at several thousand years after the formation of the first protostar. 
We should also consider the effects of turbulence and magnetic fields. In fact, the presence 
of a coherent magnetic field in primordial environment induces a strong magnetic breaking 
effect [72], leading to less fragmentation in the accretion disk. Also recent simulations predict 
more fragmentation in case with strong turbulence. Those are the issues to be addressed in 
the future. 


2.2.2. Binary black hole formation by dynamical interactions in dense star clusters. 


A. Binary black holes formed in star clusters 

Dense star clusters are one of promising formation sites of BBHs. Especially, globular 
clusters (GCs) have been expected as a source of BBHs merging within the Hubble time. 
In these days, most studies on BBH formation in GCs are being done with the Monte- 
Carlo method [7, 73]. Those works suggested a local merger rate density of ~ 5 Gpc^? yr! 
[73, 74]. They also predicted BBHs formed in GCs have several distinct features from those 
formed under isolated environments. The GC-originated BBHs have more isotropic spin- 
orbit misalignment angles [75]. Some of those BBHs may contain massive (~ 100M5) BHs 
locating in the so-called “the pair-instability mass gap" [76], or may have finite eccentricities 
at a GW frequency of — 10 Hz [77, 78]. Such Monte-Carlo simulations, however, need careful 
calibrations of some parameters from comparison with the N-body results [79, 80]. On the 
other hand, direct N-body simulations with a realistic number of particles (i.e., N ~ 10°) 
are still numerically too expensive with cost of N? and with very small binary orbital period 
compared with the GC lifetime. As a result, currently, there is only a single example of 
million-body runs performed using a GPU cluster [81], although the situation may change 
in near future by adoption of tree-direct hybrid methods [82, 83] and novel schemes for 
binary integration [84]. 

On the other hand, open clusters (OCs) have not been expected to produce such tight 
BBHs because they are less massive and less dense compared with GCs. From the mass 
function of star clusters, however, there are hundred times more OCs than GCs. Therefore, 
we thought it might be interesting to examine the formation rate of BBHs in OCs seriously. 

We performed a series of N-body simulations of OCs using a direct N-body code, 
NBODY6++GPU, which includes both stellar and binary evolution models [85]. As an initial 
condition, we adopted a Plummer model [86] with a total mass of 2.5 x 10?Mo with 0.1 solar 
metallicity. We calculated 360 realizations per a given set of the parameters and obtained 
^ 1 BBH per cluster on average. The merger time of the BBHs due to gravitational wave 
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emission is given as a function of the semi-major axis (a) and eccentricity (e) of BBHs [87]. 
In Fig. 6, the distributions of a and e of BBHs formed in the simulations are shown. We 
can see that most of BBHs distribute in red-colored region (Region 3), which correspond to 
those dynamically formed and the characteristic value of their semi-major axis is set by the 
condition that the potential depth of binaries is of the order of that of the cluster. Owing to 
shallow potential depth of OCs, their separation is too large for the BBHs to merge within 
the Hubble time. On the other hand, we found some BBHs that can merge within the Hubble 
time in blue-colored region (Region 1), in which BBHs have very low eccentricity and small 
semi-major axis. Those are formed via binary stellar evolution such as stable and unsta- 
ble mass transfer. Some of them experience a dynamical encounter after the binary stellar 
evolution and have large eccentricity (yellow-colored, Region 2). Three formation channels 
above are summarized in Fig. 6. 

In a GC, more BBHs are dynamically formed and, owing to the deeper potential well, they 
tend to have smaller value of semi-major axis and some of them merge within the Hubble 
time. On the other hand, in open clusters, the major channel for BBH mergers is different: 
a binary consisting of two massive main-sequence (MS) stars are first formed and becomes 
tight as a result of the binary stellar evolution. If this separation is kept until BBHs are 
formed after the stellar lifetime, they can merge within the Hubble time. 

From our results, we can estimate a merger rate of BBHs formed in OCs at 10 Gyr ago to 
be 0.3 yr! Gpce~? [88]. This value corresponds to 20-50% of those formed in GCs citeRo- 
driguez16. This means BBHs formed in OCs contribute to GW events much more than 
previously expected. 

We also repeat similar simulations but with different metallicities ranging from zero to one 
solar metallicity. Using the BBH formation rate and distribution of their orbital parame- 
ters, and assuming cosmic metallicity evolution with some dispersion [89] and a cosmic star 
formation rate density [90], we can estimate the merger rate of BBHs formed in OCs at 
the present universe by summing up the contributions from all OCs. The estimated event 
rate is ~ 70 yr-! Gpc ? [91], comparable to that estimated from LIGO/Virgo O2 and O3, 
64.04338 yr-1Gpc ^? [92]. 

Our simulations also give a prediction on BH-MS binary distribution observable with Gaia 
[93]. In the next data release of Gaia (Gaia DR3), binaries will be included. Using the results 
of our simulations, we estimated that ~ 10 BH-MS binaries would be detectable with Gaia 


[94], even if we assume relatively stringent observational constraints [95]. In contrast to BH- 
MS binary formed in galactic fields, those formed in OCs would have smaller mass of the 
MS companion (< 5Mo), longer orbital period (> 1 yr), and higher eccentricity (> 0.1). 


B. Binary population synthesis model of extremely low-metallicity star 

Binary population synthesis (BPS) is a powerful tool to predict merger rates of BBHs 
formed in galactic fields (e.g., [96]). Coupled with star cluster simulations, BPS can also pre- 
dict the merger rate of BBHs formed in dense stellar clusters (e.g., [74]). So far, most studies 
have focused only on BBHs formed under Pop I and II environments (0.01 S Z/Zo € 1). On 
the other hand, it has been suggested that Pop III-originated BBHs have distinct features 
from Pop I/IlI-originated BBHs [97]. Thus, it is important to bridge the metallicity gap (i.e. 
0 < Z/Zc < 0.01) in order to elucidate origins of merging BBHs. Hereafter, we call stars in 
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Fig. 6  Semi-major axis (a) and eccentricity (e) distribution of BBHs obtained from our 
simulations (from Fig. 1 of [88]). Star symbols show BBHs with a merger time less than 
10 Gyr. Colors of the symbols indicate the primary mass (Mi) of the BBHs. Colored back- 
ground indicates classification depending on their formation mechanisms illustrated in 7. 
Colored curves show the relation between a and e to give tew = 10 Gyr for BBHs with 
masses of 15, 25, 35, and 45Mo (red, yellow, green, and blue). 


this metallicity gap the extreme metal-poor (EMP) stars. In order to bridge the metallic- 
ity gap, we have made a table for evolutionary tracks of stars with —8 < log(Z/Zco) € —2. 
Note that there has been no evolution track computation for EMP stars before. We have 
compiled those evolutionary tracks as a set of fitting formulae, and incorporated them into 
BSE and NBODY6 ([98, 99], respectively). The BSE code is based on several BPS calcula- 
tions (e.g., [96, 97]), and the NBODY6 code and its derivatives are widely used for obtaining 
BBH merger rates in dense star clusters (e.g., [88, 100-105]). Here, we describe these fitting 
formulae briefly (see [106], in detail). 

We first make stellar evolution models as a reference for fitting formulae by means of 
1 dimensional (1D) simulation. For the 1D simulation, we use the HOSHI code [107-110]. 
We follow the time evolution of stars with 8 < M/Mo € 160 for —8 € log(Z/Zo) < —2. We 
start the calculation from the zero-age main-sequence (ZAMS) and stop the calculation at 
the carbon-ignition time. We do not include stellar wind mass loss. We take into account the 


17/89 


Massive MS - Massive MS Giant — Massive MS Mass Transfer BH - Massive MS 


BH - Giant Mass Transfer BH - BH 


Three-body Region 1 


Encounter co > [rJ T e *——— smaller a 
x| 


& e~0 


Mey BH - Massive MS 
Formation 
Region 3 


e > O au Ln s m larger a 


Interaction 
& Ejection 


Region 2 
smaller a 
0O<e<1 


O<e<i1 
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Fig. 8 HR diagram for comparison between stellar evolution models (black) and fitting 
formulae (red). The stellar evolution in this figure ends at the carbon-ignition time, and 
stellar wind mass loss is not taken into account. 


post carbon-ignition evolution (e.g. supernova) and stellar wind mass loss as post processing. 
This is the same as in the BSE and NBODY6 codes incorporating the evolution of Pop I/II stars. 
The evolution of EMP stars differs from that of Pop I/II in the following respects: neither 
Hertzsprung-gap (HG) nor blue-loop phase exists and the envelope of massive (M ~ 20M) 
stars remains radiative until the end of their lives. This has significant impact on binary 
evolution, such as the mass accretion rate in the Roche-lobe overflow phase, and occurrence 
of post-MS mergers and common envelope evolution. 

We use the above stellar evolution models as reference, and make fitting formulae of EMP 
stars. Figure 8 shows stellar evolution on HR diagram in the cases of log(Z/Zo) = —5 and 
—8 for comparison with the results obtained from 1D calculations and the fitting formulae. 
Note that the stellar evolution in this figure is terminated at the carbon-ignition time, and 
stellar wind mass loss is not taken into account. The tracks by fitting formulae show good 
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agreement with those by 1D calculations, with the maximum deviation less than 20%. This 
error is satisfactory for BPS and star cluster simulations coupled with BPS, considering larger 
uncertainties coming from uncertainties in tidal interaction, common envelope evolution, 
stellar wind mass loss, natal kicks of neutron star and black holes, supernova mass loss 
including pair instability, pulsational pair instability supernovae, and so on. 


2.2.3. Outlook. To distinguish the possible formation channels described above, we need 
further statistical constraints with increasing the number of BBHs identified by GW detec- 
tions. For instance, binary population synthesis models predict that the Pop III progenitor 
stars provide characteristic features in resulting chirp-mass distributions. A most prominent 
feature is a peak around 30Mo [97], which may be still consistent with the latest LIGO/Virgo 
O3a results. Future GW observations will impose more constraints on other properties of 
BBHs, such as mass ratio, separation, spin, and eccentricity, etc. Maturing both in-situ and 
ex-situ formation channels is indispensable for understanding the origins of BBHs, making 
full use of the observational constraints. 


3. Binary Neutron Stars 

3.1. Exploring NS physics 

3.1.1. BNS waveform modeling based on mumerical-relativity simulations. GWs from 
binary neutron stars (BNSs) contain rich information of the NSs. In particular, the tidal 
deformability of NSs has been proposed as one of the most promising quantities related 
to the equation of state that can be extracted from the GW observation [111-126]. The 
tidal deformability A is defined through Q;; = —A£j;j, where Qij and E;; are the induced 
quadrupole moment and the external tidal field, respectively. The tidal deformability can 
also be written in terms of the radius of the neutron star, R, and its Love number, ke, as 
A = 2k3 R9. Therefore, the behavior of \ depends on the inner structure of a NS, its EOS. 
For example, for the soft (stiff) EOSs, the radius of a NS is small (large), therefore, the NS 
is less (largely) deformed, i.e. A takes small (large) value. The simultaneous measurement 
of the mass and the tidal deformability of the NSs provides a substantial constraint on the 
equation of state of nuclear matter which is yet poorly understood [127, 149]. 

'To extract the tidal deformability of NSs from the observed signals, an accurate theoret- 
ical waveform template is crucial. For this purpose, many efforts have been made to model 
the gravitational waveforms of a NS merger in the past few decades. The waveforms includ- 
ing the linear-order tidal effects are derived by post-Newtonian (PN) calculation for the 
early inspiral stage |117, 124]. Furthermore, the waveforms employing the effective-one-body 
(EOB) formalism are derived to incorporate higher-order PN effects |111-115, 128-132]. 
However, these waveform model could not be accurate enough for the estimation of the 
tidal deformability due to a significant systematic error of the unknown higher-order PN 
terms [116, 119, 125, 126]. Furthermore, the tidal correction to the waveforms is derived 
based on the PN calculation which is only valid in the regime where the orbital separation 
is sufficiently large. Since the tidal effect is most significant just before the merger, it is 
necessary to examine the accuracy of these waveform models. 

Numerical-relativity (NR) simulation is the unique method to predict the tidal effects in 
a regime where the non-linear effect of hydrodynamics should be taken into account in the 
framework of general relativity [128, 133-142]. The recent progress in numerical techniques 
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and computational resources enabled us to perform the high-precision NR simulations for 
BNS mergers and to generate the waveforms with systematic errors sufficiently small for the 
GW data analysis. 

In our recent work |140, 143, 144], we perform NR simulations for NS-NS mergers system- 
atically varying the mass and the tidal deformability of NSs: 6 and 4 models that cover the 
symmetric mass ratio, 7, in the range of 0.247—0.25, the binary with the total mass of ~ 2.7 
and 2.5 Mo, and 5 different EOSs which cover the binary tidal deformability in the range of 
e: 300-1900. Here, the binary tidal deformability, A, is defined by 


A= Š [( + 7-317?) (Ar + Aa) 


— V1 - 4n (14-95 — 117”) (Ay — I] ; (3) 


where A; 2 A;/m? (i — 1,2) are the tidal deformability of each individual star (Note that 
mi < m»). By this work the waveforms from NS-NS mergers for more than 15 inspiral orbits 
are obtained for a wide range of binary parameters. 

We carefully estimate the error budgets of the waveforms obtained by the NR simulations. 
In our NR simulation code, the field equations are solved under discretization. We found 
that the error due to the finite grid resolution is the largest error budget among the error 
sources of the numerical simulation. To estimate the error due to the finite grid resolution, we 
perform the simulations with various gird resolutions and checked the convergence property 
of the obtained waveform data. As the results, we find that the phase error of the waveforms 
is a sub-radian order [140, 143, 144]. 

Employing the obtained high-precision NR. waveforms, we examined the validity of the 
analytical waveform models for NS binary. We show that the latest tidal-EOB (SEOB- 
NRv2T [129, 132]) waveforms can be accurate even up to ~% 3ms before the onset of 
merger [140] (see Fig. 9 for an example). However, the phase difference between the tidal- 
EOB waveforms and the NR results is still larger than ~ 1rad after two NSs come into 
contact for the case that the NS radii are larger than ~% 13 km. This finding indicated that 
further improvement of the waveform model would be needed to suppress the systematic 
error in the measurement of the tidal deformability. 

For this purpose, we develop a phenomenological model for GWs from inspiraling NSs 
based on our high-precision NR. simulations. Our phenomenological waveform model is 
derived in the frequency domain as in the Phenom-series for binary black holes [145] for 
convenience in data analysis. We decompose a frequency-domain gravitational waveform, 
h (f), into the frequency-domain amplitude, A (f), and phase, V (f), (with an ambiguity in 
the origin of the phase) by 


h(f) - AQ)e 70, (4) 
and we define the corrections of the NS tidal deformation to GW amplitude and phase by 
Atiaal (f) = A (f) — Assn (f) (5) 
and 
Vaga (f) = Y (f) — Yesu (F), (6) 
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Fig. 9 Comparison between NR and SEOBNRv2T waveforms for the models with the 
total mass of ~ 2.5 Mc and the symmetric mass ratio of 0.247. The top and bottom panels 
show the amplitude of NR and SEOBNRv2T waveforms and the phase difference between 
them, respectively. Here, D and mo denote the distance to the source and the total mass of 
the binary. The vertical dashed lines are the peak amplitude time when the GW amplitude 
reaches a peak in each model. 


respectively. Here, Aggy (f) and Uppy (f) are the GW amplitude and phase of a binary 
black hole with the same mass as the BNS, respectively (hereafter referred to as the point- 
particle parts). In this work, we employ the SEOBNRv2 waveforms [146] as the fiducial 
point-particle part of GWs. 

For the tidal part phase and amplitude of our phenomenological waveform model, we 
employ the following function forms motivated by the 2.5 PN order post-Newtonian formula 
Refs. [115] and [139] for the tidal-part amplitude and phase correction: 


tidal _ 3 | 39; (1 i TDI "P 


~ 1289 | 2 
3115 3/2, 28024205 , 4283 s/z 
1 ! 7 
j ( * 1248" 3302208 ^ — 1092^" (7) 


for the phase correction and 


5r] me ~ _ 
Atidal [EU "0 A 7/4 
V24 Da 


27 5 4496 1. 
8 
«( ig? T ba) (8) 


for the amplitude correction. Here, a, p, b, and q are the free parameters of the models. 
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Fig. 10 (Left) Difference in the tidal-part phase between the hybrid waveforms and our 
analytical model for the various binary configurations (see [140] and [143] for the details 
of the models). Phase differences are plotted after the alignment in the frequency range of 
10-1000 Hz. (Right) Relative difference of amplitude between the hybrid waveforms and our 


analytical model. 


Since our NR waveforms only contain the waveforms of which frequency is larger than 
x 400 Hz, they are too short to calibrate the phenomenological waveform model only by 
themselves. Thus, for this purpose, we construct hybrid waveforms employing our NR wave- 
forms and the effective-one-body waveforms of Refs. [129, 132] (SEOBNRv2T) as the high 
and low-frequency parts, respectively, and use them to calibrate the phenomenological wave- 
form model. The hybridization of the waveforms is performed in the time-domain by the 
procedure described in Ref. [139, 143]. 

'To focus on the inspiral-phase waveform and to avoid the contamination from the post- 
merger waveforms, which would have large uncertainties, we restrict the GW frequency range 
in 10-1000 Hz. The fitting parameters are determined by employing the hybrid waveforms 
with the largest value of binary tidal deformability in the models studied in this work. By 
performing the least square fit with respect to the phase difference and relative difference 
of the amplitude, we obtained a = 12.55, p = 4.240, b = 4251, and q = 7.890. Employing 
these parameters, we find that the phase and amplitude corrections to the frequency domain 
waveforms are reproduced within 0.1 rad and 25% for all the binary configurations studied 
in this work (see Fig. 10). 

To validate our phenomenological waveform model more quantitatively from the view point 
of the data analysis, we calculate the mismatch between the our waveform models and the 
hybrid waveforms, F', defined by 


( A | js etri Hee] 
[lea] {Pall 


F = 1 — max 
Qo.to 


where (-|-) and ||- || are defined by 
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Fig. 11 ‘Tidal phase in the frequency domain normalized by the leading, Newtonian 
(relative 5PN-order) tidal phase formula. Our analytical model (referred to as Kyoto tidal) 
assuming A = 1000 (dot-dashed, blue), 400 (dashed, blue), and 100 (dotted, blue), the 
NRTidal model (solid, red [135]) and the NRTidalv2 model (solid, cyan [136]) are shown in 
the plot. Note that latter two models are independent of A when normalized by the leading 
tidal phase. 


and 
p= Vll = f (Afà). (11) 


Here Sa denotes the one-sided noise spectrum density of the detector, and we employ 
the noise spectrum density of the ZERO DETUNED HIGH POWER configuration of advanced 
LIGO [147] for it. We found that the mismatch between our phenomenological waveform 
model and the hybrid waveforms is always within 10^? [143, 144]. 

Gravitational waveform models for BNSs based on NR waveforms are also derived in [135] 
and [136] in a similar manner. The main difference between our and their work is the dif- 
ference of the NR waveforms and the tidal-EOB waveforms used for the model calibration. 
Figure 11 compares our tidal phase correction model and those of [135] and [136] normalized 
by the leading-order PN term (see also [148] and [136]). Figure 11 shows that those models 
agree with each other within 10-20% for 10-1000 Hz. Indeed, those models give consistent 
results for the current detector sensitivity as shown in the next subsection (see also [136]). 


3.1.2. Measuring tidal deformability from GW. As mentioned above, NSs provide unique 
laboratories to study the properties of ultra-dense matter. To constrain the equation of 
state of NSs effectively, we need accurate waveform models for the data analysis. We assume 
that the waveform of the GWs can be decomposed into the point-particle, spin and tidal 
parts. The point-particle and the spin parts are the same as those of the binary black hole 
case, on the other hand, tidal part is that of the BNS case because of the matter effects 
of NSs. The standard waveform, so-called TaylorF2 (hereafter TF2), is derived by the post- 
Newtonian (PN) approximation [152]. The point-particle [153, 154], and spin parts [155-157] 
are calculated up to the 3.5PN-order, while the effect of the tidal part (let us call it PNTidal) 
enters from the 5PN-order, calculated up to the 7.5 PN order [115]. Since the tidal effects 
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become important at the late inspiral stage where the PN approximation becomes invalid, 
we need better waveform models to estimate tidal deformability accurately. 

Dietrich et al. [135, 136] (NRTidal) and Kawaguchi et al. [143] (KyotoTidal) independently 
construct new waveform models in which the tidal part is improved by the calibration using 
NR waveforms. In KyotoTidal, nonlinear effects of A is taken into account. The lack of the 
point-particle and spin parts beyond the 4PN-order also bias the measurement of binary 
quantities. Kawaguchi et al. have also constructed higher PN-order correction terms of the 
point-particle part by fitting the SEOBNRv2 waveform model [146, 158], so-called TF2+. Here, 
TF2+ KyotoTidal is calibrated in the frequency range of 10-1000 Hz, because they focus 
on the inspiral phase to avoid uncertainties in the post-merger phase. Several studies have 
also derived constraints on the NS EOS via measuring tidal deformability from GW170817 
[159-169]. 

Narikawa et al. separately analyze the GW data of a binary-neutron-star merger, 
GW170817, by each of the Advanced LIGO twin detectors [170]. They found that the poste- 
rior probability distributions of the binary tidal deformability for the Hanford and Livingston 
detectors are distinctively different as shown in Fig. 12. They have suggested that the dif- 
ference between the detectors cannot be understood by statistical error only, since only the 
posterior with the Livingston detector does not change smoothly as the upper cutoff fre- 
quency increases. Their findings suggest that further research of the noise properties in the 
high-frequency region of the Livingston data is needed to extract information on the tidal 
deformability of GW170817. Narikawa et al. also found that there is a difference in the esti- 
mates of A for GW170817 between NR calibrated waveform models (the KyotoTidal and 
NRTidalv2 models) as shown in Fig. 13 [148]. Here, NRTidalv2 model is an upgrade of the 
NRTidal model [136]. The order of peak values of A for the different waveform models in 
Fig. 13 can be understood by the difference in the magnitude of the tidal phase at around 
A e 600 — 900 shown in Fig. 11. However, the difference is smaller than the statistical error. 

The second BNS merger event GW190425 was observed during O3a [171]. The total mass 
of the system was about 3.4Mo, which is larger than the total mass of any other known BNS 
system. Since the LIGO Hanford detector was offline at the time of the event, parameter 
estimation was performed by using the data from the LIGO Livingston and Virgo detectors. 
The EOS with large tidal deformability A > 1100 was disfavored, which is consistent with 
the several analyses done for GW170817. For the low-spin prior (component spin is enforced 
as x13 ~ U[—0.05, 0.05]), the 90% upper limit of GW190425 becomes much smaller than 
that of GW170817, A < 600. Since the total mass of the system was large, the possibility of 
BH-NS system is discussed in [171]. 


3.2. Association of a short-duration gamma-ray burst 


3.2.1. High-energy Observations. The observations of GW170817 by high-energy instru- 
ments, Fermi, MAXI, Swift, CALET and Chandra, are reviewed in this subsection. 


Fermi. On August 17, 2017, at 12:41:06.5 (UTC), hard X-ray emission from GRB 
1708174 triggered Fermi Gamma-ray Burst Monitor (GBM with a coverage of 8keV—40MeV; 
[172]) near the detection limit [173]. The emission with a short hard pulse and a soft tail 
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Fig. 12  Marginalized posterior probability distribution of binary tidal deformability, A 
for GW170817, derived by data of different detectors with fmax = 2048 Hz. (Reprinted with 
permission from [170]. (C)(2017) by the American Physical Society.) 
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Fig. 13  Marginalized posterior PDFs of binary tidal deformability, A for GW170817, 
estimated by different waveform models for fmax = 1000 Hz. 


lasted for —2 s, as shown in Fig. 14. From the obtained time duration, this GRB is classi- 
fied as a short GRB. Also the SPectrometer on board INTEGRAL Anti-Coincidence Shield 
(SPI-ACS) detected a hard pulse completely coincident with the Fermi/GBM one [174]. 
Amazingly, the binary coalescence detected by the Laser Interferometer Gravitational- 
wave Observatory (LIGO) occurred ~1.7s before the Fermi/GBM trigger [9] and the sky 
position determined by the LIGO GW observation is consistent with the Fermi/GBM one. 
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Fig. 14 Light curve of GRB 170817A detected by GBM in the 8-50keV and 50-300keV 
bands. Emission consists of a hard pulse and a weak soft tail. The vertical solid line repre- 
sents the time of the binary neutron star merger. The horizontal dashed line represents the 
averaged background level. 


The gamma-ray spectrum for the initial hard pulse is well fitted by a power-law function 


with a photon index of a = -1.60.4 and an exponential cutoff at Epeak = 185 + 62keV. For 
the soft tail emission, a blackbody with kgT = 10.3 + 1.5keV well represents the spectrum. 
The fluence over the total duration is ~3 x 107 "erg cm" ?, which falls in the ~50th percentile 
of the fluence distribution obtained by Fermi/GBM [173]. Assuming a distance to the host 
galaxy NGC 4993 as 43Mpc, the isotropic equivalent energy and peak luminosity in the 
lkeV — 10MeV band are ~3 x 10*°erg and ~2 x 10*’ergs~', respectively, which are by 3-6 
orders of magnitude smaller than those of other short GRBs detected by Fermi/GBM. The 
observed features for GRB 1708174 would suggest that the hard X-ray emission came from 


a misaligned jet to the observer (i.e., off-axis jet). 

The higher-energy gamma-ray observation of GRB 170817A by Fermi Large Area Tele- 
scope (LAT; [175]) covering the 0.1-300GeV band was unfortunately not available due to 
entry to the South Atlantic Anomaly at the time of occurrence of the merger time [17]. 
Fermi/LAT resumed scientific operation ~10°s after the merger time and put a flux upper 
limit of 4.5 x 10^ Pergcm ? s^! (95% confidence level) in the 0.1-1GeV band using the 
initial 1000s observation after the operation resumption. 

Motivated by study of GRB 1708174, investigation of previous short GRBs with a simi- 
lar signature has been performed [176-179]. In particular, the nearby short GRB 150101B 
detected by Fermi/GBM (z = 0.134 [180]), which could accompany the kilonova and the 
off-axis jet emission suggested by the optical and X-ray observations [177], has a hard spike 
and a soft tail. Thus, the two-component sigunature could be a common feature for short 
GRBs [176]. The isotropic equivalent energy and luminosity in the 1keV — 10MeV band 
are ~2 x 10“erg and ~4 x 10°°ergs~!, respectively, which are in a fiducial range of short 
GRBs. One of the possible interpretation for these phenomena suggests that GRB 150101B 
could originate from a more on-axis jet than GRB 1708174. 
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MAXI. Monitor of All-sky X-ray Image (MAXT, [181]) is an instrument for monitoring 
X-ray sky, which was launched in 2009 and has been observing from the International Space 
Station (ISS). Because MAXI does not have a moving mechanism and is mounted to the 
ISS, we cannot actively control the pointing direction. The orbital period of the ISS is 92 
minutes and thus an X-ray source is usually observed once in the same period. Since MAXI 
covers about 85% of the sky every orbit, the source activity can trace back to the time before 
the event. 

In studying X-ray counterparts of GW events, we use the Gas Slit Camera (GSC; [182]) 
GSC consists of 12 position sensitive proportional counters. Six of them covers 3 degree x 160 
degree field of view (FOV) toward the horizontal direction and the others compose zenithal 
FOV with the same dimension. These FOVs move as the ISS rotates. A scan observation 
(transit) for a point source lasts 40-150 seconds depending on the source position in the 
FOV [183]. 

The effective area of the GSC for a source changes with time linearly like a triangle. For a 
source with constant flux (or that can be assumed as constant in a scan transit), we calculate 
an average flux by dividing the source count by the effective exposure, that is a sum of the 
effective area along the time. 

For GW170817 MA XI did not detect the X-ray emission. The upper limits are summarized 
in [184]. Fig. 15 shows the luminosity and upper limits for early (< 10 days) X-ray afterglows 
[185, 186] in the case of GW170817. The position of the EM counterpart lay at the gap 
between the horizontal and the zenithal cameras for the first three orbits and the first scan 
started after 4.7 hours. Moreover the first three observations were incomplete and the upper 
limits were 1-2 order of magnitude higher than the usual. However, it is important to mention 
that this was the earliest observation of the X-ray counterpart of the GW event. 

Although the observation start time was late in this case, if MA XI's observation starts ~ 
100 seconds after the trigger, MAXI could observe an X-ray afterglow equivalent to the level 
observed by Swift at several hours after the trigger. 

'The best lesson of the observation of GW170817 was that it was important to continuously 
observe the sky as long as possible. In O3 we have increased the observation efficiency to 
reduce the gap in the observation. 


Swift. The Neil Gehrels Swift Observatory (Swift; [187]) consists of three scientific instru- 
ments. The Burst Alert Telescope (BAT; [188]) is a wide field of view telescope monitoring 
about 1/6 of the sky in the 15-150keV band. The X-ray Telescope (XRT; [189]) is an X-ray 
telescope covering the energy range in the 0.3-10keV band in 24' field of view. The UV- 
Optical Telescope (UVOT; [190]) observes the wavelength range of 170-600nm in 17' field 
of view. 

At the merger time of GW170817, the error region of GW170817 was occulted by the 
Earth and not visible by BAT [185, 191]. Therefore, no useful information can be derived 
from the BAT data. The XRT and UVOT started their observations about an hour after 
the merger time around the center of the Fermi GBM error region (~ 1? radius). After the 
localization information became available from the LIGO and Virgo at 4.8 hours after the 
merger time, Swift started a series of 120s observations at the positions of known galaxies 
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Fig. 15 The observed luminosity and upper limits of the X-ray counterpart of GW170817 
in early time. MAXI upper limits are based on the result of [184]. The first six points are 
for each scan transit and the last two are for one day and 10-day observations. The data of 
Swift XRT and Chandra are converted to the luminosity at 2-10keV assuming the spectral 
parameters in their original papers [185, 186]. The dotted, dashed, and dash-dotted lines 
are models for luminosity L as a function of the time t, L c t^, and a are -1, -1.5, and 
-2, respectively. The horizontal gray line is a typical GSC upper limit for one scan transit, 
scaled to the distance of GW170817 (40Mpc). 


inside the 33.6 deg? error region of the gravitational-wave detectors. No new X-ray sources 
(2107 erg cm? s^!) were found during those observations [185]. 

Swift started the observations at the position of the optical transient of GW170817 about 
14.4 hours after the merger time. Although the XRT found no X-ray from the position of 
the transient, the UVOT detected a fast decaying UV emission which is interpreted as the 
blue kilonova associated with GW170817. The non-detection of X-ray at the early phase by 
XRT places a crucial difference to a typical short GRB afterglow which is dominated by an 
on-axis emission [185]. 


CALET. The Calorimetric Electron Telescope (CALET; [192, 193]) is the scientific 
instrument attached on the International Space Station (ISS). The CALET Gamma-ray 
Burst Monitor (CGBM;[194]) is the gamma-ray burst monitor onboard CALET and is 
capable to observe emissions from 7keV to 1MeV (Hard X-ray Monitor; HXM, FoV of 
^60? from the boresight) and from 40keV to 20MeV (Soft Gamma-ray Monitor; SGM, 
FoV of —110? from the boresight) utilizing two different detectors. No significant emission 
was observed by CGBM at the time of the merger of GW170817. The 9096 upper limit is 
1.3 x 10~“erg cm~? s7! 
ISS [36]. 

In addition, the main CALET instrument, the high-energy Calorimeter (CAL), observes 


in the 10-1000keV band assuming no shielding by the structure of 


gamma-rays from ~1GeV up to 10TeV with a field of view of nearly 2sr. The field of view 
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Fig. 16 The X-ray light curve of the counterpart of GW170817 observed by Chandra. 
'The source count-rates are converted to the unabsorbed flux in the 0.3-8keV band using the 
Galactic absorption column Ng = 7.6 x 102? cm"? and the photon index of 1.59 [198]. 


of CAL did not include the location of GW170817 at the time of the merger. A search 
for delayed emission over the 2-month period following the event resulted in 9096 C. L. 
upper limits on the energy flux of 1.2 x 107 !Perg cm? s^! for gamma-rays above 1GeV and 
4.0 x 107 Perg cm"? s^! above 10GeV [195]. 


Chandra. The Chandra X-ray Observatory, which provides the best sensitivity in X-ray 
for a point source, first targeted to the optical transient of GW170817 at 2.2 days after the 
merger time. Although no X-ray emission was detected in this first visit, Chandra detected 
the significant X-ray emission on the 2nd visit at 8.9 days after the merger in the isotropic 
luminosity of 9 x 10?9ergs^! at a distance of 40Mpc [196]. 

According to the year-long monitoring observations by Chandra, the X-ray flux continues 
to increase up to ~160 days after the merger time by a power-law index of 0.9 [197, 198]. 
After the peak, the flux shows the decline by a power-law index of ~ —2 (Fig. 16). There 
is a hint of an X-ray flaring feature around 155 days after the merger time [199]. Chandra 
has been collecting the data up to 2.5 years after the merger and continues monitoring. This 
unique X-ray light curve behavior along with the radio data will be discussed in the following 
section. 


3.2.2. Theoretical Interpretation. The GW event GW170817 [9] and the associated short 
GRB 1708174 [36, 173, 174] have provided, for the first time, clear evidence for a relativistic 
jet from a binary neutron star merger observed from off-axis (as schematically shown in 
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Fig. 17 Electromagnetic counterparts associated with the off-axis jet from short GRB 
170817A (from Fig. 1 of Ref. [13]). 


Fig. 17). The observations also show that the jet is structured, excluding a top-hat distri- 
bution of energy and Lorentz factor of the jet. We also discuss the current issues about the 
exact structure of the jet, the origin of the gamma-ray emission, and future prospects. 


A. Evidence of an off-axis jet 

The afterglow observations of GW170817A verify the existence of a relativistic jet: the 
superluminal motion of the radio image indicates a relativistic collimated source [11, 12] and 
the closure relation between the spectral index and the light curve slope after the luminosity 
peak is consistent with the standard afterglow model of a relativistic jet [198, 200, 201], 
excluding a failed-jet scenario that the jet is choked by the ejecta from the neutron star 
merger [202-204]. 

The jet should be strong enough for successfully breaking out the merger ejecta in this 
event [205, 206]. The breakout time since the jet launch (t = to) is approximately estimated 
as 


1/2 -1/2 =i 
to — tm Liso 0 Liso 0 
ty — to ~ 0.17 i 0.078s | ——— 12 
PIN 5 ( 0.1s ) (= wo) x 5 (z ergs-! U2) 


where Liso,o is the isotropic luminosity at the base of the jet and to — tm is the time delay 
of the jet launch since the merger (modifying Eq. (32) of Ref. [207]). Here the second term 
is necessary for this event because the ejecta expansion velocity is comparable to the jet 
head velocity, in contrast to collapsars with a static envelope. In order to have the short 
GRB 1708174 ~ 1.7s after the merger, the jet luminosity should be large enough Liso,o Z 
3 x 10Pergs^! and the delay time should be short enough to — tm € 1.3 s [207]. 

'The observed isotropic energy of the gamma-ray emission is very small, several orders 
of magnitude smaller than typical [36, 173, 174]. Therefore the jet should be off-axis for a 
reasonable radiative efficiency [13, 14]. An off-axis observer receives photons emitted outside 
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Fig. 18 Left: Examples of the jet structure that are consistent with the afterglow obser- 
vations of GRB 1708174: a hollow-cone (magenta line) and Gaussian jet (dashed line). The 
green region is not in principle determined by the current data. Right: The corresponding 
light curves for a hollow-cone (thick dashed lines) and Gaussian jet (thin dashed line) with 
the data. The blue line is used for the inverse reconstruction (modifying Fig. 8 of Ref. [217]). 


the relativistic beaming cone, and the apparent energy of the off-axis jet becomes faint (see 
also below). 


D. Diverse jet structures 

The afterglow observations, in particular the slowly rising light curves from radio to X- 
ray, are not consistent with a top-hat jet [208], but indicate a structured jet [12, 198, 201, 
209—214]. These are the first confirmation of previous suspicion that realistic GRB jets are 
structured [215, 216]. 

However, there remains confusion about the jet structure as different authors give different 
structures, such as Gaussian [198, 201, 211] and power-law profiles [12, 209] (see Fig. 1 of 
Ref. [13]). These structures are obtained by assuming a functional form of the jet structure 
with model parameters, and adjusting the parameters through fitting the theoretical light 
curves to the observations. 

We instead invent a novel method to determine the functional form of the jet structure itself 
[217]. This method solves an inverse problem, uniquely reconstructing a jet structure from 
a given off-axis GRB afterglow by integrating an ordinary differerential equation, which 
is formulated based on the standard theory of GRB afterglows. Applying the inversion 
method to GRB 1708174, we clarify that the current observational errors still allow various 
jet structures, discovering by-product that a hollow-cone jet is also consistent with the 
observations as shown in Fig. 18. 

An important caveat is that the current afterglow observations only constrain the jet core 
(magenta in Fig. 18), not in principle the outer jet structure, (green; which is crucial for 
the gamma-ray emission as discussed below). Early afterglow observations are required to 
obtain the outer jet structure [217]. 

'The origin of the jet structure is still unclear. First, the jet may have intrinsic structure 
from the beginning of the launch. Second, the jet may be structured during the propagation 
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Fig. 19 The surface brightness distribution of the prompt emission from a Gaussian jet 
(from Fig. 3 of Ref. [14]) 


if the outer part loads baryon from the ejecta. Third, the cocoon component (the collided 
jet and ejecta) may contribute to the outer structure since a part of the cocoon is relativistic 
due to partial mixing of the baryon. Early macronova/kilonova observations are important 
for probing the cocoon, which has crucial information on the jet. 


C. Origin of the gamma-ray emission 

The prompt gamma-ray emission of GRB 1708174 generally comes from an off-center jet, 
neither the jet core nor the line of sight but the middle, as in Fig. 19 [14], given the jet 
structure with a rapidly decaying tail as in Fig. 18. This is because the emission from the jet 
core is suppressed by the relativistic de-beaming, while the line-of-sight emission becomes 
too faint to dominate the off-axis emission from the inner jet. The off-center emission solves 
seemingly puzzles that the observed prompt spectrum is inconsistent with the spectral Amati 
relation [218, 219] (because the off-center is different from the usually-observed jet core) and 
causes the compactness problem [220, 221] (because the off-centre jet is much less energetic 
and much closer to the line of sight than the jet core). 

'The off-center structure still has large uncertainties because it is not constrained by the 
observed afterglows (as shown by the green region in Fig. 18). Whether it is a jet or cocoon is 
also unknown. The observed gamma-rays might be the jet emission scattered by the cocoon 
[222]. Future observations of off-axis events will reveal the jet structure, where roughly ~ 10% 
events are expected to be brighter at smaller viewing angles than GRB 1708174. 


3.8.  Optical/Infrared Counterparts and Heavy Element Nucleosynthesis 

'The origin of the elements in the Universe is a subject of interest in astronomy and astro- 
physics. In particular, the origin of the elements synthesized by rapid neutron capture, or 
so-called r-process, is one of the long-standing, unsolved problems. Theoretically, NS merg- 
ers are one of the promising candidate sites of r-process nucleosynthesis [223-231]. However, 
there has been no observational evidence of r-process by NS mergers. 
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Fig. 20 Optical, infrared, and radio telescopes in the J-GEM network. 


When NS mergers occur, a small fraction of mass of NSs is ejected into interstellar space 
[228, 230, 232-235]. In the ejected material, r-process nucleosynthesis is expected to take 
place. Then, we expect optical and infrared emission powered by radioactive decays of freshly 
synthesized nuclei [236—241]. This phenomenon is called “kilonova” or *macronova" (here- 
after we call it kilonova). In other words, if we can detect kilonova after detection of GWs 
from a NS merger, we can test r-process nucleosynthesis by the NS merger. 

The event rate of NS mergers can be measured from GW observations. The ejection of 
r-process elements per event can be estimated from electromagnetic (optical/infrared) obser- 
vations. Therefore, by these multi-messenger observations, we can study the production rate 
of r-process elements by NS mergers and test if NS mergers can be the origin of r-process 
elements in the Universe. 

'To achieve observations of kilonova, we coordinated the Japanese collaboration for Gravita- 
tional wave ElectroMagnetic follow-up (J-GEM, Section 3.3.1). Thanks to this observational 
network, we succeeded in observations of the electromagnetic counterpart of GW170817 
(Section 3.3.2). By combining our state-of-art numerical simulations (Section 3.3.3), we 
advance our knowledge about the physics of NS mergers and the origin of r-process elements 
in the Universe. 


3.3.1. J-GEM. J-GEM is a network of optical, infrared, and radio telescopes for EM 
follow-up observations of GW sources. Figure 20 shows a summary of telescopes in the 
J-GEM. The roles of the optical telescopes can be divided into two categories: one is the wide- 
field surveys (Kiso, MOA-II, and Subaru Hyper Suprime-Cam (HSC)) and galaxy targeted 
surveys (other telescopes). For more details of each telescope, see Morokuma et al. (2016) 
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Fig. 21 Left: Localization map of GW170817 and Subaru/HSC pointing [265]. Right: 
three-color image of the counterpart of GW170817 [267]. 


[242]. J-GEM signed Memorandum of Understanding with LIGO/Virgo collaboration in 
2014, and started to receive the alert of GW detection from the first observing run of 
Advanced LIGO (O1). 

We demonstrated coordinated observations using this observing network for the first GW 
source GW150914 (BBH, [1]). Within 4 days from the GW detection (2 days from the initial 
alert), about 24deg? area has been observed with Kiso Schmidt telescopes. Also, within 
6-7 days from the GW detection, 18 nearby galaxies were observed with the B&C 61-cm 
telescope. Although no EM counterpart has been identified, this was the first attempt of 
multi-messenger observations of GW sources [242]. 

More intensive observing campaign has been performed for the second GW source 
GW151226 (BBH, [243]). For this event, 238 nearby galaxies were observed by galaxy 
targeted surveys [244], which is great improvement compared with the observations of 
GW150914. Also, Kiso Schmidt telescope and MOA-II telescopes covered 778deg? and 
145deg? area, respectively [244]. Furthermore, Subaru/HSC covered 63.5deg? area with 
unprecedented sensitivity (limiting magnitudes of 24.6 and 23.8 for the i and z bands, 
respectively, [245]). The total area covered by the wide-field surveys was 986.5deg?, which 
corresponds to ~29% of the probability map of GW151226. It is worth emphasizing that 
the sensitivity of Subaru/HSC observations was deep enough to detect expected kilonova 
emission even at about 200Mpc, which is about the design sensitivity of Advanced LIGO, 
Advanced Virgo, and KAGRA. 

Subaru/HSC is one of the most powerful wide-field imager: its field-of-view (1.8deg?) is 
the largest among 8-10m class telescopes. Table 3 summarized our follow-up observations 
using Subaru/HSC. Details of follow-up observations of $190510g are presented by Ohgami 
et al. (2021) [277]. 
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Table 3 Subaru/HSC observations for GW sources 

! Detected with low-letency pipeline but not recovered with the offline analysis. 
2 Initially reported as BNS with a high probability but can be terrestrial noise. 
8 Initially reported as Mas Gap object but later classified as BBH. 


Target Type Disc. date Obs. date Filter Area Probability Lim. mag 
(UUT) (UT) (deg?) (76) 

GW151220 BBH 2015-12-26 2016-01-07 HSC-i 63.1 12 24.4 
2016-01-07  HSC-z 63.3 12 24.0 
2016-01-13 HSC-i2 63.1 12 24.4 
2016-01-13 HSC-z 63.2 12 24.0 
2016-02-06 HSC-i2 65.0 12 24.5 
2016-02-06 HSC-z 65.0 12 24.0 
G275404.  BH-NS 2017-02-25 2017-02-27 HSC-z 26.0 0.5 23.2 
GW170814 BBH 2017-08-14 2017-08-16 HSC-Y 32.5 35 22.2 
2017-08-17 HSC-Y 32.6 35 22.2 
GW170817 BNS 2017-08-17 2017-08-18 HSC-z 262 61 21.7 
2017-08-19 HSC-z 27.3 61 21.6 
$190510g BNS? 2019-05-10 2019-05-10 HSC-Y 118.8 1.1 22.9 
$191216ap BBH? 2019-12-16 2019-12-20 HSC-z 2.3 0.035 25.2 
8$200224ca BBH 2020-02-24 2020-02-25 HSC-r2 57.8 80.8 25.3 
HSC-z 57.8 80.8 23.6 
2020-02-28 HSC-r2 57.8 80.8 24.6 
HSC-z 57.8 80.8 23.2 
2020-03-23 HSC-r2 57.8 80.8 25.5 
HSC-z 57.8 80.8 23.8 


3.9.28. GW170817. The first GW detection from a NS merger has been achieved for 
GW170817 [9]. The detection triggered EM follow-up observations covering entire wavelength 
range [10]. Thanks to the observations with three GW detectors, the position of the GW 
source is localized to ~ 30deg?. In optical/infrared wavelengths, a counterpart AT2017gfo 
was identified by several groups [246-268]. 

We have performed wide-field survey using Subaru/HSC, covering 23.6deg? [265]. The 
observations recovered the detection of AT2017gfo. Furthermore, thanks to the wide coverage 
and good sensitivity, we statistically ruled out that other transients in the localization area 
are associated with GW 170817, or in other words, we concluded that AT2017gfo is the most 
likely counterpart. 

Intensive optical and infrared observations of AT2017gfo have been performed in the 
framework of J-GEM network [267]. Observations with Subaru/HSC show that the z-band 
brightness quickly declines. On the other hand, the observations with IRSF/SIRIUS shows 
that near-infrared light curves evolve more slowly. This is fully consistent with the expected 
behavior of kilonova. 

Using radiative transfer simulations, we find that the observed light curves can be explained 
by ^ 0.03Mo of ejecta including lanthanide elements, which have a high opacity [269]. Also, 
the blue optical component, revealed by observations with Subaru/HSC, B&C/Tripol5, and 
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Fig. 22 Density distribution and electron fraction distribution in the post-merger ejecta, 
obtained by viscous-radiation hydrodynamics simulations [271]. 
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Fig. 23  Multi-color light curves of GW170817/AT2017gfo (circles) compared with the 
results of two-dimensional radiative transfer simulations (solid lines) [272, 274] 


MOA-II/MOA-cam3, requires an ejecta component with a smaller fraction of lanthanide. 
Since the estimated total ejecta mass is higher than the expected ejection by the first dynami- 
cal ejection from the NS merger, we concluded that the observed signals are mainly produced 
by post-merger mass ejection. 
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3.3.8. Advances of Theoretical Models. The observed properties of AT2017gfo are con- 
sistent with the expectation from numerical relativity simulations [270]. In particular, to 
have an ejecta component with a moderate lanthanide fraction, at least temporal presence 
of hypermassive neutron star is suggested. If the merger results in the prompt collapse to 
a BH, there is no source of neutrino radiation, and thus, almost entire ejecta are expected 
to be lanthanide-rich, which is not consistent with the observations. In this way, numerical 
relativity simulations play crucial roles to connect the initial conditions estimated from GWs 
(mass and mass ratio) with the final outcome of kilonova. 

Mass ejection from NS mergers can be roughly divided into two phases: (1) dynamical 
mass ejection mainly by the tidal force and shock interaction in the dynamical time scale 
and (2) subsequent post-merger mass ejection from the accretion torus. Observations of 
GW170817/AT2017gfo highlighted the importance of post-merger mass ejection from NS 
mergers. However, the detailed simulations have been difficult. We have performed long-term 
general relativistic neutrino radiation hydrodynamics simulations, by taking into account the 
effects of viscosity [271]. We find that post-merger ejecta (> 107 ?M) dominate over the 
dynamical ejecta (< 10 27M). Thanks to the inclusion of neutrino transfer, we could also 
provide the distribution of electron fraction in the ejecta. We showed that the ejecta have 
high electron fraction (> 0.25), which means that post-merger ejecta is lanthanide-poor 
(Fig. 22). 

Depending on the mass and mass ratios of the NS merger, various ejecta components 
with different properties (mass, velocity, and electron fractions) are expected. Therefore, it 
is important to study the observable outcome using the results of numerical relativity. We 
have performed multi-dimensional radiative transfer simulations by consistently taking the 
interplay of multiple ejecta components into account [272]. We showed that AT2017gfo can 
be naturally explained by the geometry predicted by numerical relativity simulations (Figure 
22). 

We then extended our radiative transfer simulations for various initial conditions. Using 
the updated atomic data [273], Kawaguchi et al. (2020) [274] showed that the properties 
of kilonova should have diversity depending mainly on the total mass of the system. Since 
this is multi-dimensional simulations, we can also predict the variety of the luminosity for 
different viewing angles. It is shown that optical emission is greatly suppressed when NS 
merger is observed from the equatorial plane while infrared emission is almost unchanged. 
These simulations will be useful to facilitate future EM follow-up observations. 


3.3.4. Future Prospects. Our observational and theoretical results indicate that NS merg- 
ers synthesize a wide range of r-process elements. With the current estimate of the event 
rate and mass ejection, NS merger can explain the total amount of r-process elements in our 
Galaxy [3, 4]. However, the event rate is still largely uncertain. Also, we still do not fully 
understand that the mass ejection is universal or not. Varieties in the mass ejection and 
nucleosynthesis are expected by numerical simulations. Therefore, it is important to observe 
more NS mergers to accurately derive the event rate and to obtain the general picture of 
the mass ejection and nucleosynthesis. 

The sensitivity of GW detectors are still improving. In the third observing run (O3), a 
typical sensitivity corresponds to the range of the NS merger up to 130Mpc (LIGO Liv- 
ingstone), 110Mpc (LIGO Hanford), and 45Mpc (Virgo). Thanks to these high sensitivity, 
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Fig. 24 Expected number of NS merger as a function of distance assuming the event rate 
estimated from the O2 run [5]. The dashed lines show the range for 9096 confidence level. 


two reliable GW events including at least one NS have been reported. However, due to the 
poor positional localization or large distances, any EM counterpart has not been identified. 
We refer the reader to Sasada et al. (2021) [278] and Ohgami et al. (2021) [277] for our 
observing effort in O3. The situation will be improved when the sensitivity of Virgo becomes 
comparable. Furthermore, when KAGRA joins the network with a comparable sensitivity, 
positional localization will be greatly improved (down to a few deg? for the signals such as 
GW150914, [275]). 

Figure 24 shows the expected event number per year as a function of the distance. Once all 
the sensitivity of GW detectors reach ~ 200Mpc, we expect > 10 events every year. Kilonova 
associated with such events is expected to be faint, ~ 22 mag in optical within 2-3 days after 
the merger. Therefore, observations with Subaru/HSC will play important roles. Also, deep 
near-infrared imaging and spectroscopic observations with large-aperture telescopes, such as 
Subaru, Keck, Gemini, and upcoming TAO 6.5m telescope [276], will be important for firm 
identification of kilonova. Such multi-wavelength observations will provide us with unique 
information of r-process nucleosynthesis. 


4. Supernovae 

4.1. GW signatures from core-collapse supernovae 

4.1.1. The CCSN GW signatures. Unlike the GW signals from compact binary coales- 
cence where a template-based search is best suited (e.g., [279]), gravitational waveforms 
from CCSNe are essentially of stochastic nature. T'his is because the waveforms are affected 
by turbulence in the postbounce core, which is governed by the multi-D, non-linear hydro- 
dynamics. In order to clarify the GW emission mechanisms, extensive numerical simulations 
have been done so far in different contexts (e.g., [280-287] and [288] for a review). For canon- 
ical supernova progenitors [289], core rotation is generally too slow to affect the dynamics 
(e.g., [290, 291]). For such progenitors, the GW emission takes place in the postbounce phase, 


which is characterized by prompt convection, neutrino-driven convection, proto-neutron star 
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(PNS) convection, the standing accretion shock instability (SASI, [292]), and the g(/ f)-mode 
oscillation of the PNS (e.g., [287, 293-297]). 

The most generic GW emission feature seen in recent self-consistent 3D CCSN models 
is the one from the PNS oscillation [298-301]. The characteristic GW frequency increases 
almost monotonically with time due to an accumulating accretion to the PNS, which ranges 
approximately from ~ 100 to 1000Hz. On the other hand, the typical frequency of the SASI- 
induced GW signals is concentrated in the lower frequency range of ~ 100 to 260Hz and 
persists when the SASI dominates over neutrino-driven convection [299-302]. The detection 
of these distinct GW features is considered as the key to infer which one is working more 
predominantly in the preexplosion supernova core, namely neutrino-driven convection or the 


SASI [299]. 
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Fig. 25 An example gravitational waveform of the + (blue line) and x (orange line) mode 
from a 3D-GR CCSN model of 15Mo star (model SFHx in [300]). The time is measured after 
core bounce (:t = 0). A source distance of 10kpc is assumed. The waveform was extracted 
via a standard quadrupole formula with GR corrections [285]. Since the waveform is weakly 
dependent on the source orientation, an unbiased direction is chosen (e.g., from the north 
pole ((0, 9) = (0,0)). This figure is taken from [304]. 


Shown Fig. 25 is an example waveform taken from a 3D CCSN simulation in [300]. For 
the model, the hydrodynamics evolution is self-consistently followed in full general relativity 
(GR), starting from the onset of core-collapse of a 15Mo star [303], through core bounce, 
up to ~ 350ms after bounce. As consistent with the outcomes from recent 3D models (e.g., 
(299, 301]), the hydrodynamic evolution and the associated GW waveform is characterized 
by the prompt convection phase shortly after bounce (Tyy S; 20ms with Tpp the postbounce 
time, shown simply as "time" in Fig. 25), then the linear (or quiescent) phase, which is 
followed by the non-linear phase (7, Z 140ms) when the vigorous SASI activity (as well as 
the growing GW amplitudes) was observed for the model (see [300] for more details). The 
dominance of the SASI over neutrino-driven convection persists over 140ms < Typ S; 300ms, 
after which neutrino-driven convection dominates over the SASI. 

In order to discuss the detectability of the CCSN GW signals (like of Fig. 25), previous 
studies have traditionally relied on a GW spectrogram analysis, aiming to specify the GW 
feature by the excess power in the time-frequency domain (by taking the square norm of the 
short-time Fourier transform). However, the spectrogram analysis has a trade-off relation 
between the time and frequency resolution, leading to a limited time-frequency localiza- 
tion. Because of this drawback, some features could have been potentially overlooked in the 
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previous spectrogram analysis in the context of CCSN GWs. Recent developments in the 
time-frequency analysis (TFA) have yielded various time-frequency representation alterna- 
tives to the spectrogram. In particular, quadratic time-frequency representations, such as 
the Wigner-Ville distribution and its modified forms, enable us to perform high-resolution 
TFA (e.g., [305-307]). 

Recently, the S-method utilizing the Wigner- Ville distribution was applied in [304], which 
successfully extracted several key (CCSN) GW features in the time-frequency domain (left 
panel of Fig. 26, where the waveform (Fig. 25) was used). In comparison with the conven- 
tional TFA based on the Fourier transform (e.g., Fig. 1 in [300]), the S-method significantly 
improves the sharpness of the modes. Furthermore by defining the instantaneous frequency 
(IF) of the excess in the time-frequency domain (corresponding to the bright-color regions 
in the left panel of Fig. 26, see [304] for more details), five modes (A B, C, C*, and D, the 
right panel of Fig. 26) were identified. 

Mode A in the right panel of Fig. 26 corresponds to the PNS g(/f)-mode oscillation, which 
is excited in the vicinity of the PNS surface by the downflows and/or directly originates from 
the deceleration of infalling convective plums [297]. Modes B and C are quasi-static modes, 
in the sense that these two modes barely change with time in the time-frequency domain. 
Mode B (fg ~ 130Hz) is originated from SASI, the frequency of which is twice of the SASI 
frequency (fsAsr ~ 65Hz) because of the quadrupole nature of the GW emission (see also 
[302]). Mode D is the overtone of the mode B (fp © 2fp), which is most likely to come from 
the PNS core oscillation (see Fig. 8 in [304]). Mode C intersects mode A at Ty ~ 180ms, 
after which this mode is denoted as C*. Although yet to be investigated elaborately, modes 
C and C* are likely to arise in the innermost region (~ 10km) [308]. Note that the typical 
frequency of the SASI-induced GWs is in the range of ~100 — 260Hz (e.g., modes B and D 
in the left panel of Fig. 26), which is in the best sensitivity range of the currently running 
interferometers. For a Galactic event, the signal-to-noise ratio (SNR) of the GW signals 
predicted in the most recent 3D CCSN models was estimated in the range of — 4 — 10 
for the advanced LIGO [299]. With the third-generation detectors on line (e.g., Cosmic 
Explorer and Einstein Telescope), these signals would be never missed for the CCSN events 
throughout the Milky way. 


4.1.2. Circular Polarization of CCSN GWs. Recently, yet another GW feature has been 
reported, which is circular polarization of the CCSN GWs. The importance of detecting the 
GW circular polarization was first pointed out by [309] in the context of rapidly rotating 
core-collapse. More recently, Hayama et al. (2018) [310] presented analysis of the GW circular 
polarization using results from core-collapse of a non-rotating 15 Mo star [300], the waveform 
of which is already shown in Fig. 25. 

Figure 27 visualizes the time evolution of the GW circular polarization (hereafter, CP) of 
models SFHx (left panel) and TM1 (right panel) from [310]. Note that SFHx [311] and TM1 
denotes the name of the nuclear equation of state (EOS) employed in the 3D-GR simulation 
[300]. T'he key difference is that the softer EOS (SFHx) makes the PNS radius and the shock 
radius at the shock-stall more compact than those of TM1. This leads to much stronger 
activity of the SASI for SFHx compared to TM1. 

In fact, one can see a clearer polarization signature for SFHx (left panel of Fig. 27) charac- 
terized by the bigger GW amplitude with the right-handed (red line) and left-handed mode 
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Fig. 26 The left panel is the time-frequency representation of the sample waveform 
(Fig. 25), which is extracted by the S-method. The right panel shows the GW mode identi- 
fication of the spectrogram based on the instantaneous frequency identification method via 
the S-method (see [304] for more details). In the right panel, the white dashed line repre- 
sents the theoretical prediction of the peak GW frequency of the PNS g-mode oscillation 
[297]. Note that the deviation of mode A from the white dashed line (especially in the later 
postbounce phase) was previously observed in the literature [287, 299]. These figures are 
taken from [304]. 


=29 -22 
5 x10 . 5 x10 . 
di 1f 
= 0 0 E y 
-1 L =4 L 
E : : : B ; 
-2 E 0 1 2 -2 E O 1 2 
h, x10777 h, x10 


Fig. 27 Shown are the trajectory of the GW circular polarization on the hy — hx plane 
of models SFHx (left panel) and TM1 (right panel) of 3D-GR simulations using different 
equations of state (see text). The two characteristic epoch with the right-handed (140ms € 
Tpb S; 200ms) and left-handed (220ms S$ Tpb S; 320ms) polarization are highlighted with the 
red and blue color. T'he gray line denotes the whole trajectory during the simulation time. 
The source distance is assumed at 10kpc. These figures are taken from [310]. 


(blue line) than those for TM1 (right panel). Note in each panel that the two SASI-dominant 
phases of SFHx are colored by red or blue (see also bottom panels of Fig. 1 in [310]). The 
non-axisymmetric flows associated with spiral SASI give rise to the circular polarization. 
Before Tpb ~ 140ms, the CP is small because the SASI activity is still weak. But, after 
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Tpb ~ 140ms when the (spiral) SASI activity begins to be vigorous (e.g., [300]), the right- 
handed CP firstly emerged, which is followed by the left-handed CP (the blue line) until 
Ty ~ 320ms, after which neutrino-driven convection dominates over the SASI. 

To explore the detectability of the CP, Monte Carlo simulations of the coherent net- 
work analysis [312] were performed for the signal reconstruction, where the network of 
LIGO Hanford, LIGO Livingston, VIRGO and KAGRA was considered. They reported the 
enhancement of the signal-to-noise ratio of the GW circular polarization (SNRcp) relative 
to that of the GW signal itself (SNRrr). Although this is likely because the Gaussian noise 
has little component of the circular polarization, more detailed analysis is required to draw 
a robust conclusion whether or not the GW CP could provide a new probe to clarifying the 
CCSN inner-workings such as the SASI and the PNS oscillations (see [310] for more detail). 


4.1.8. Rapid Rotation: GW and Neutrino signals from an exploding 21? Mc, star. Rapid 
rotation in the iron core (with the initial rotation frequency typically greater than 0.5rad/s) 
leads to significant rotational flattening of the collapsing and bouncing core, which produces 
a time-dependent quadrupole (or higher) GW emission (e.g., [313] for a review). For the 
bounce signals having a strong and characteristic signature, the iron core must rotate enough 
rapidly!. The GW frequency associated with the rapidly rotating collapse and bounce is in 
the range of ~ 600 — 1000Hz [280]. A current estimate based on a coherent network analysis 
using predictions from a set of 3D models shows that these GW signals could be detectable 
up to about ~ 20kpc for the most rapidly rotating model ([312], see also [315, 316]). 

Figure 28 shows the neutrino and GW signals [317] obtained in 3D core-collapse supernova 
simulation of a rapidly rotating 27Mo star that is exploding due to the growth of the so- 
called low-T/|W | instability [290]. The time modulation seen in the left panel corresponds to 
the neutrino light-house effect where the spinning of strong neutrino emission regions around 
the rotational axis leads to quasi-periodic modulation in the neutrino signal. Depending on 
the observer's viewing angle, the time modulation will be clearly detectable in IceCube and 
the future Hyper-Kamiokande. The GW emission is also anisotropic where the GW signal is 
emitted, as previously identified (see [288] for a review), most strongly toward the equator at 
rotating core-collapse and bounce, and the non-axisymmetric instabilities in the postbounce 
phase lead to stronger GW emission toward the spin axis. The right panel in Fig. 28 shows 
that these GW signals can be a target of LIGO-class detectors for a Galactic event. The origin 
of the postbounce GW emission (e.g., the bar-mode (m — 2) deformation of the PNS due 
to the low T/|W| instability), naturally explains the reason that the peak GW frequency 
is about twice of the neutrino modulation frequency. These results demonstrate that the 
simultaneous detection of the rotation-induced neutrino and GW signatures could provide 
a smoking-gun signature of a rapidly rotating PNS at the birth (see also [320, 321]). 


4.1.4. Hydrodynamics and GW signals from a BH-forming massive star. Successful obser- 
vations of GWs from binary black holes (BBHs) by the LIGO-VIRGO collaborations are now 
posing a new challenge; what is the origin of the massive BHs ? One of the most plausible 


! Although rapid rotation and the strong magnetic fields in the core are attracting great attention 
as the key to solve the dynamics of collapsars and magnetars, one should keep in mind that recent 
stellar evolution calculations predict that such an extreme condition can be realized only in a special 
case [314] (SS 1% of massive star population). 
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Fig. 28 The left panel shows detection rates of v, at 10kpc for a rapidly rotating 27 Mo 
model as a function of time after bounce [317]. The red and green lines correspond to the 
event rates (per 1ms bin) for an observer along the equator and the pole (e.g., parallel to the 
rotation axis), respectively. A quasi-periodic modulation (red line) corresponds to ~ 120Hz 
(red lines) seen for the observer along the equator. The right panel shows the characteristic 
GW spectra relative to the sensitivity curves of advanced LIGO, advanced Virgo ([318]) 
and KAGRA ([319]). The peak GW frequency (the vertical red band) is about twice of the 
neutrino modulation frequency. 


scenarios is a binary stellar evolution in a low-metallicity environment (see [6] for a review). 
It has been proposed that two massive stars in the approximate range of 40 to 100Mo lead 
to the formation of a massive helium core after experiencing the Roche-lobe overflow and 
common envelope phase (e.g., [322—324] for collective references therein). The gravitational 
collapse of the massive core (~ 30M) could account for some of the relevant BH mass 
ranges (at least in the high-mass end) in the GW events, although the formation path to 
the massive core and further to the BH is still very uncertain due to the complexity of the 
binary evolution and the fallback dynamics (e.g., [325]). 

In order to clarify the formation process of the BH, one requires GR neutrino radiation- 
hydrodynamics core-collapse simulations of such massive stars in 3D. Due to the high 
numerical cost, most of the previous studies with BH formation have been done assum- 
ing spherical symmetry (1D) (e.g., [25, 27] and collective references therein). In the context 
of multi-D simulations with multi-energy neutrino transport, Ref. [326] reported 1D and two- 
dimensional (2D) core-collapse simulations of a solar-metallicity 40Mọ star using two-flavor 
IDSA scheme ([327]) and a post-Newtonian gravity to include GR effects. More recently, a 
BH-forming 3D-GR simulation of a zero-metallicity 40M. star with an approximate neu- 
trino transport (FMT) scheme was reported in [328]. It is only recently that the first 3D-GR 
simulation (using a 70Mọ star [329]) with detailed neutrino transport [330] following the 
dynamics up to BH formation was published [331]. In the simulation, the evolution equations 
of metric, hydrodynamics are solved based on the Baumgarte-Shibata-Shapiro-Nakamura 
formalism [332, 333] and the evolution of neutrino radiation field based on the multi-energy 
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Fig. 29 The left panel shows time evolution of the maximum rest-mass density (solid 
line) and minimum lapse function (dotted line) for a 70Mọ star (denoted as Z70, red line) 
and a 40Mo star (as S40, black line), respectively. The right panel shows angle-averaged 
profiles of the (rest-mass) density (black line) and the ratio rs/R (red lines, see text) for Z70 
as a function of the enclosed baryon mass for three representative time slices near the final 
simulation time. "-0.321" and ”-2.359” ms denote the time before the final simulation time 
(Tän = 293ms after bounce). The green line corresponds to rs/R for S40. These figures are 
from [331]. 


M1 scheme [331]. Casting the helium core of ~ 31Mo, the zero-metallicity 70M. star was 
chosen as one of the possible progenitor candidates of the BBHs. 

The left panel of Fig. 29 shows temporal evolution of the maximum density pmax (solid 
lines) and the minimum lapse Qin (dotted lines) for the 70M star (denoted as Z70, red 
line) and the 40M, star (as S40, black line), respectively. Note that S40 is shown for compar- 
ing with previous results. The compactness parameter of S40 is much smaller (£2.5 ~ 0.26) 
than that of Z70. The maximum density at bounce (pmax ~ 4 x 10!^g cm7?) is quite similar 
between Z70 and S40. After bounce, the increase of the maximum density of Z70 (red solid 
line) is significantly faster than that of S40 (black solid line). In both Z70 and S40, the min- 
imum lapse (dotted lines) shows a gradual decrease after bounce. At around Tg, — 293ms 
after bounce, it shows a drastic drop to Qmin = 0.0645 for Z70, indicating that the PNS core 
starts to collapse rapidly toward a BH formation. 

The right panel of Fig. 29 explains how 75, is related to the BH formation time. Shown in 
the panel is the profiles of (angle-averaged) density (black lines) and a diagnostics to measure 
the BH formation as a function of the enclosed baryon mass My at some representative 
snapshots near Tg, for Z70 (red lines) and at Tpb = 400ms for $40 (green line). As a BH- 
formation diagnostics, the ratio of r;/ R is shown, where rs and R denotes the Schwarzschild 
radius and the radial coordinate, respectively. One can see that the maximum r,/ R is ~ 0.3 
at 2.359ms before Tg, (the thin red line labelled by ”-2.359 [ms|") and rapidly increases with 
time, approaching to unity (precisely, 0.932) at Tg, (thickest red line), which was judged 
as the epoch of the BH formation in [331]. It should be noted that for the unambiguous 
definition of the BH formation, one requires the implementation of the so-called apparent 
horizon finder (e.g., [334]) in numerical relativity simulation, which is still left as a future 
work. 

At the (fiducial) BH formation time, the mass and the radius is My.) BH ~ 2.60(2.51) Mo 
and Riso ~ 4km, respectively. By contrast, S40 shows significantly less compact structure 
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Fig. 30 A snapshot of the entropy distribution (kg baryon !) for Z70 at Tp, ~ 294ms just 
before the BH formation ([331]). The sheet represents the lapse function (a) on the z — 0 
plane. This figure is taken from [331]. 


(green line) at the final simulation time (Tpp = 400ms). The BH formation should occur much 
later, possibly when the mass shell at R(My = 2.6Mo) ~ 10?cm accretes to the stalled shock. 
Using the same EOS (LS220 [335]), this expectation is in line with [326], which reported the 
BH formation at T,, ~ 700ms and with [328] at Ty ~ 1 s. 

The time evolution of the (angle-averaged) shock radii Rs, the gain radius Rg, the ratio of 
the advection timescale to the neutrino-heating timescale in the gain region Tagy/Theat, and 
the mass in the gain region Mgain are presented in Fig. 2 of [331], respectively. For model Z70, 
the shock revival was obtained after Tjj Z 260ms. At this time, the maximum temperature 
in the core becomes as high as T ~ 100 MeV at a slightly off-center region at Riso ~ 10km 
(equivalently at My ~ 1.0Mọo). Subsequently the high temperature region propagates out- 
ward in the mass coordinate, although spatially inward, due to the continuous mass accretion. 
The maximum temperature reaches ~ 170 MeV at Rigo ~ 1km (My ~ 1.4Mọ) just before 
Tán. In this second collapse phase to the forming BH, the high neutrino emission makes 
the heating timescale shorter than the competing advection timescale in the gain region. 
Aided by strong convection behind the shock, the stalled shock is revived at Ty, Z 260ms 
(Tadv/Theat 2 1). This also results in the increase in the gain mass (see the blue line) due to 
the shock expansion. 

Figure 30 visualizes a 3D hydrodynamics feature near at the BH formation. During the 
first ~ 160ms after bounce, the neutrino heating is still weak and high entropy bubbles are 
yet to appear. Only after Tpb Z 230ms, the formation of high entropy plumes (s = 15kg 
baryon-!) was seen due to the intense neutrino heating from the hot PNS. At this time, 
the mass in the gain region Mgain also starts to increase. The expansion of the (merging) 
high entropy plumes was observed, leading to the shock revival. The lapse function shows 
a steepest drop in the center (see the cusp in the plane of Fig. 30), which corresponds to 
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Fig. 31 A gravitational waveform for the BH-forming core-collapse of a 70M, star. Note 
that h} and D denote the GW amplitude of the + polarization and the distance to the 
source, respectively. 


the BH formation. By expanding the shock radius into the spherical harmonics, we find that 
the deviation of the shock from spherical symmetry (in the low-modes / — 1,2) is less than 
~ 2%. This clearly indicates that neutrino-driven convection dominates over the SASI in 
this case. 

Finally, Fig. 31 shows the GW prediction for the Z70 star. The waveform is extracted 
along positive z-axis via the quadrupole formula. The GW amplitude (multiplied by the 
distance to the source) stays at small value < 10cm before the shock expansion occurs 
(Tpb ~ 260ms). The strong GW emission thereafter mainly originates from strong convection 
motion behind the shock. Just before the final simulation time, the GW amplitude reaches ~ 
100cm, though only for a short duration. Although more quantitative discussion is apparently 
needed, a rough estimate shows that the GW signal could be targeted by third-generation 
GW detectors such as the Einstein Telescope (ET) and the Cosmic Explorer (CE) ([433, 434]) 
for the source distance less than ^ Mpc. Regarding the neutrino signals, both electron 
type (anti-)neutrinos show decreasing and plateau trend for Tjj Z 260ms, whereas heavy- 
lepton neutrinos show a rapid increase both in the luminosity and energy. These features 
are consistent with those previously identified in 1D full-GR simulations with Boltzmann 
neutrino transport ([336]), which is due to rapid contraction of the PNS to the forming BH 
(see also, [337, 338] ). The detection of the short-live (~ 300ms after bounce) neutrino signals 
are basically limited to Galactic events (see [32] for a review). However, further study would 
be needed to clarify the contribution of these BH forming massive stars to the prediction of 
diffuse neutrino supernova background (e.g., [339, 340]). 


4.1.5. Short Summary and Future Prospects. To summarize, each phase of a CCSN has a 
range of characteristic GW signatures and the GW polarization that can provide diagnostic 
constraints on the evolution and physical parameters of a CCSN and on the explosion hydro- 
dynamics. The GW signatures described in this section commonly appear in the frequency 
range of 100 - 1000Hz in recent simulations of the CCSN. Among them, the PNS oscillatory 
modes are currently recognized as the model-independent GW signature. The characteristic 
frequency (fp x M/R?) of the fundamental (f/g) modes is predominantly dependent on the 
PNS mass (M) and the radius (R) [341]. In order to break the degeneracy, the detection 
of other eigen-modes (p, w modes) of the PNS oscillations is mandatory. In fact, the GW 
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asteroseismology of the PNS has just started [342-344], the outcomes of which should reveal 
the requirement of the next-generation detectors to detect the whole of the eigen-modes 
(presumably extending up to several kHz) for the future CCSN event. 

Regarding the SASI-induced GWs for a Galactic supernova source, they are predicted to 
be detectable by the advanced LIGO with the singal-to-noise ratio (SNR) in the range of 
~ 4-10 by the most recent 3D CCSN models [299]. With the third-generation detectors (e.g. 
Cosmic Explorer and Einstein Telescope) online, these signals would never be missed for 
CCSN events throughout the Milky Way. A current estimate of SNR for GWs from rapidly 
rotating collapse and bounce based on predictions from a set of 3D models using a coherent 
network analysis shows that these GW signals could be detectable up to about ~ 20kpc for 
the most rapidly rotating model ([312]; see also Refs. [315, 316]). See a recent review by Ref. 
[345] for a more comprehensive review regarding the detectability of the CCSN GW signals 
by the future detectors. 

After decades of progress, the CCSN theory is finally converging, but not without detailed 
numerical effort and much theoretical sweat. But at last not in the least, we have to draw a 
caution that the current numerical results and the associated GW predictions that we have 
reported in this Section should depend on the next-generation calculations by which more 
sophistication can be made not only in the neutrino transport schemes (such as Boltzmann 
transport [346—348]), but also in the treatment of GR with the BH formation. Therefore 
we provide here only a snapshot of the moving (long-run) documentary film that records 
our endeavours for making our dream of " GW astronomy of the next Galactic CCSN” come 
true. 


4.2. Understanding Supernovae via Neutrino Emissions 


One of the primary sources of near-field GWs, core collapse supernova explosions are among 
the most dramatic and important events to take place in the Universe. Agents of great 
destruction — anything within tens of light-years is utterly annihilated by the dying star — 
they are nevertheless responsible for the existence of life itself, for they (and the collisions of 
neutron stars they produce) are the sole source of all elements heavier than helium. Therefore, 
understanding these complex explosive events is necessary if we are to understand why we 
are here, and why the Universe looks the way it does today. 

Supernova neutrinos, famously observed from SN1987A, provide a unique and vital probe 
into the inner dynamics of these events. Released together with GWs during the initial stel- 
lar collapse, neutrinos and GWs are both certain to travel through any obscuring dust or 
gas and remain undiminished upon their arrival at the Earth. Neutrinos also carry infor- 
mation regarding the end state of the star: for explosions within our galaxy, collapses into 
neutron stars or black holes, the eventual sources of far-field GWs, can be differentiated via 
observations of neutrino emissions. 

Our goal is to make theory and experiment work together so that we will be ready to make 
the best possible observations of the next galactic explosion and maximize our extraction of 
information on the explosion mechanism, progenitor, and nuclear physics ingredients. 

In order to both predict and make sense of the neutrino signals from the next explosion 
in the Milky Way galaxy, the world's most advanced supernova numerical simulations are 
needed. Using nuclear data such as equation of state and neutrino reactions, and by solving 
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the 6D Boltzmann equation in 2D/3D, we can provide detailed predictions of the time 
distribution and energy spectra of supernova neutrinos. 

To verify these complex computational models, we will enrich the famous Super- 
Kamiokande (SK) detector with gadolinium [Gd] salt. Doing so will turn it into the world’s 
most advanced supernova neutrino detector, capable of real-time tagging and identification 
of individual supernova neutrino interactions with nanosecond-scale time resolution. This 
will make the diffuse supernova neutrino flux from all past supernova explosions visible for 
the first time. Starting around the middle of 2020, a gadolinium-loaded SK will collect a 
steady stream of supernova neutrino data, the first such new data in over 30 years. 

The new experimental measurements made possible by gadolinium will then be fed into 
the theoretical models, testing their predictions using real supernova neutrino data and 
allowing us to better prepare for the next nearby supernova explosion. The models will be 
refined as needed and as indicated by the past supernova data. The gadolinium-loading of 
SK will also greatly improve the detector’s response to a nearby supernova. Therefore, both 
experimentally and theoretically, we will be well-prepared and ready for it. 


4.29.1. SK-Gd. Water Cherenkov detectors such as Kamiokande, IMB, and SK, have been 
used for decades as effective detectors for neutrino interactions and nucleon decay searches. 
While many important measurements have been made with these very large detectors, a 
major drawback has been their inability to efficiently detect the presence of thermal neutrons. 

If a water Cherenkov detector could be improved to observe the neutrons produced by the 
inverse beta process (Zp — etn), then backgrounds would be greatly reduced. As a result, 
supernova neutrinos from explosions 35,000 times more distant than SN1987A could be seen 
by an improved SK, covering 5096 of the entire Universe. Perhaps more important for less 
distant explosions, neutron tagging of inverse beta events would facilitate the de-convolution 
of a galactic supernova's various signals, allowing a much more complete interpretation of 
the physics of the burst. In coincidence with a GW signal, detailed information concerning 
the explosion's dynamics would become even more valuable. 

The key is to add 0.2% by mass of a soluble gadolinium compound like gadolinium sulfate, 
Gd2(SOx4)3, to the water. Doing so will make >90% of the neutrons visible as a conse- 
quence of the gamma rays released by gadolinium’s capture of thermalized neutrons. This 
technique and the various new scientific advances it would make possible was first pro- 
posed in the Physical Review Letters article, “GADZOOKS! Anti-neutrino spectroscopy 
with large water Cherenkov detectors"[356]. The publication of this paper introduced the 
concept of gadolinium-enhanced water Cherenkov detectors, a transformational technology 
with a strong impact on the physics community. It provided a clear and cost-effective path 
to extend and build upon the pioneering Kamioka work of years past. 

In order to test this new technology, a dedicated gadolinium R&D facility was built in 2009 
in the Kamioka mine near SK. Called EGADS (Evaluating Gadolinium's Action on Detector 
Systems) [357], it consists of a gadolinium-loaded 200-ton scale model of SK, complete with 
240 50-cm photomultiplier tubes, its own DAQ and readout electronics, a novel selective 
water filtration system, and water transparency evaluation equipment. Based on its successful 
operation, in 2015 the SK Collaboration formally approved the plan to add gadolinium to 
SK, designating this new phase of the experiment “SK-Gd”. The T2K Collaboration, a 
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long-baseline neutrino oscillation experiment based at J-PARC which uses SK as their far 
detector, formally approved the gadolinium-loading plan in 2016. 

In order to prepare SK for the addition of gadolinium, the detector had to be opened and 
refurbished for the first time in twelve years. There were four main tasks: 


(1) Replace the photomultiplier tubes that had failed (a few hundred out of 13,000) since 
the previous in-tank refurbishment in 2006. 

(2) Fix a small water leak in the SK tank. 

(3) Clean any rust and other dirt that had accumulated in the detector since its original 
completion in 1996. 

(4) Install additional water piping to increase the total water flow for increased water 
purification, and to enable better control of the flow direction in the tank. 


This in-tank work, which required over 3000 person-days of effort, was successfully carried 
out between May 2018 and January 2019. By February 2019 the detector had been refilled 
with pure water, and was ready for the addition of Gd2(SO4)3. After taking the T2K beam 
schedule into account, the first gadolinium is now expected to go into the tank in mid-2020. 
By searching for the small but constant, diffuse flux of neutrinos produced by all core collapse 
explosions since the onset of stellar formation in the early Universe, it is expected that by 
2022 we will have collected the world’s first additional supernova neutrinos since SN1987A. 


4.2.9. Theoretical studies of supernova explosion and neutrino emission. In the rest of this 
section, we report the progress of theoretical study of supernova neutrinos for the detection of 
SK with the emphasis on detailed description of neutrino transfer and microphysics. In order 
to provide the prediction of supernova neutrino detections at SK with gadolinium loading, 
it is mandatory to get rid of the uncertainties in numerical simulations to describe the 
explosion mechanism and neutrino emission. One of the remaining uncertainties of supernova 
simulations is the neutrino transfer, which has been routinely treated with approximate 
ways in multi-dimensional simulations. We perform the sophisticated numerical simulations 
of neutrino-radiation hydrodynamics based on the direct solver of Boltzmann equation. The 
first principle-type calculations enable us to determine the final outcome of explosion, provide 
the solid prediction of neutrino emissions and evaluate the uncertainties of microphysics such 
as equation of state. 


4.29.8. Boltzmann-radiation-hydrodynamics simulation of the core-collapse supernovae. 
We have performed several simulations using the Boltzmann-radiation-hydrodynamics code 
under 2D/3D. The details of the code are illustrated in [358-360]. Various nuclear equations 
of states (EOSs), progenitor models, and the rotational velocities are employed. Here, we 
show the important features of the obtained results. 

First, we examined the effect of the nuclear EOSs [361]. We employ the Lattimer-Swesty 
(LS) EOS [362] and the Furusawa-Shen (FS) EOS [363]. The former (latter) EOS model 
employs the soft (hard) nuclear force model and the single nuclear approximation (nuclear 
statistical equilibrium treatment) for the nuclear composition. The progenitor model is the 
11.2 Mo progenitor model taken from [364]. The simulations with the different EOS models 
show a significant difference: the LS model shows shock revival while the FS model not 
(the left panel of Fig. 32). A useful quantity to understand the difference is the timescale 
ratio, which is defined as Taqy/Theat; Where Tady :— Mgain/ M and Theat := Egain/Qgain are the 
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Fig. 32 (Left) The time evolution of the shock radii for the LS (red) and FS (blue) models. 
'The thick solid lines are the angular averaged shock radii, and the dashed lines show the 
ranges between the minimum and the maximum shock radii. (Right) The time evolution of 
the timescale ratio for the LS (red) and FS (blue) models. 


advection timescale and the heating timescale, respectively. The advection timescale is the 
timescale with which a fluid element flows through the gain region. The heating timescale is 
the timescale with which a fluid element is heated and become gravitationally unbound. The 
gain region is the region where the neutrino heating exceeds cooling. The symbols Mgain, M, 
Egain, and Qgain are the mass in the gain region, the mass accretion rate, the total energy 
including the gravitational binding energy in the gain region, and the net neutrino heating 
rate in the gain region, respectively. If the timescale ratio exceeds unity, the neutrino heating 
proceeds sufficiently fast and the explosion succeeds. The right panel of Fig. 32 shows that 
the timescale ratio of the LS model exceeds unity, while that of the FS model not. Although 
M, Egain, ANd Qeain are similar between the LS and FS models, Mgain for the LS model 
is larger than that of the FS model. This is because the turbulence is stronger for the 
LS model. Indeed, the prompt convection is stronger for the LS model, and it enhances 
the neutrino driven convection at the later stage. The stronger prompt convection in the LS 
model is originated from the stronger photodissociation of the accreting nuclei: the accretion 
flow outside the shock in the LS model contains more heavy nuclei, and hence the energy 
consumed by the photodissociation of these nuclei is larger for the LS model. T'herefore the 
proper treatment of the nuclear composition of an EOS is important to assess the influence 
of EOSs on the CCSNe. 

Second, we investigated the effect of rotation [365]. The employed EOS and progenitor 
model are the FS EOS and the 11.2 Mọ model from [364]. The rotational velocity profile 
is so-called the shelluler rotation profile: Q(r) = 1rad/s/(1 + (r/1000 km)?). However, the 
employed rotational velocity is too slow to affect the postbounce dynamics even though the 
velocity is almost the highest according to the current stellar evolution theory [366]: the 
shock radii and other dynamical features are similar between the rotating and non-rotating 
models. What is more interesting here is the momentum space distributions of neutrinos. 
Figure 33 shows the neutrino angular distributions in the laboratory frame outside the 
shock on the equator. Especially, the distribution of the high-energy neutrinos is tilted to 
the ¢-direction. This is because matter rotates and drags the neutrino distribution to the 
rotational direction. This detailed angular distribution is only accessible by the Boltzmann 
solver. Since we have such information, we can assess the accuracy of the M1-closure scheme, 


50/89 


e; 


1MeV 
4MeV 
19 MeV 


€p 


Fig. 33 The angular distributions of neutrinos outside the shock on the equator at 12 ms 
after bounce. The left panel shows the 3D distributions, and the right panel shows the 


sections by the r-@ plane. The different colors correspond to the different energies as indicated 
in the left panel. 
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Fig. 34 The eigenvalues of the Eddington tensors calculated from the distribution function 
(solid) and the M1-closure scheme (dashed) in different directions at 12 ms after bounce. The 
largest eigenvalue is called the Eddington factor (red), and the other eigenvalues are named 
lateral 1 and 2 (blue and green). The upper row shows the eigenvalues itself, while the lower 
row displays the fractional difference between the solid and dashed lines. 


one of the approximate neutrino transport method. The second angular moment of the 
distribution divided by the zero-th moment is called the Eddington tensor. The M1-closure 
scheme estimates the Eddington tensor from the energy flux and density of neutrinos. In 
Fig. 34, we compare the eigenvalues of the Eddington tensor calculated from the distribution 
function directly and the M1-closure scheme. The largest eigenvalue, or the Eddington factor, 
is ~ 20% larger for the M1-closure scheme. This difference originates from so-called ray- 
collision between outward and inward rays. If we can get some information about these rays 
from neighboring matter, we may improve the accuracy of the M1-closure scheme. 
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Third, we suggested a new mechanism for the proto-neutron star kick motion [367]. We 
performed the core-collapse simulation with the 15 Mo model of [364] and Togashi-Furusawa 
(TF) EOS [368]. The TF EOS is one of the most realistic EOSs. The most interesting feature 
of this model is the proto-neutron star kick. Thanks to the exact treatment of the proto- 
neutron star in our code, the proto-neutron star moves with the velocity of O(10) km/s. So 
far, the driving force of the proto-neutron star motion is considered to be the gravity of the 
asymmetric ejecta [369]. However, the driving force observed in our simulation seems the 
recoil of the asymmetric neutrino emission. With the proto-neutron star motion, the asym- 
metric distribution of neutrinos are sustained, and hence the driving force by the neutrinos 
persists. T'his result suggests the different mechanism of the proto-neutron star motion. 

Finally, we have performed the 3D simulation with the  Boltzmann-radiation- 
hydrodynamics code. Due to the limited computational resources, the dynamics and the 
neutrino distributions until ~ 20 ms after the core bounce are investigated. With this simu- 
lation, the neutrino distribution function with no spatial symmetry is obtained for the first 
time in the world. This result provides us the important clues to understand the behavior 
of the neutrino in the supernovae. 


4.2.4. Light curve of supernova neutrinos. In order to extract the information as much 
as possible from the detection of supernova neutrinos from the next supernova in future, 
it is necessary to prepare the templates of neutrino signals by systematically covering the 
variation of progenitors and microphysics. One of such systematic sets is the supernova 
neutrino database provided by [370]. This data is now routinely used to evaluate the event 
rates at SK, replacing the rather classic data set by [371]. 

'The supernova neutrino database is constructed by the combination of supernova sim- 
ulations for the set of progenitors with different masses and metallicities. The numerical 
code of general relativistic neutrino-radiation hydrodynamics under the spherical symme- 
try [372, 373], which is the first-principle calculation, is utilized to study the series of 
phenomena; gravitational collapse, core bounce, formation of central object, the shock prop- 
agation and the associated neutrino emission. The numerical code for the thermal evolution 
of proto-neutron star with the flux-limited neutrino diffusion in general relativity [374, 375] 
is utilized to study the neutrino emission from the cooling phase of the proto-neutron star 
over 20s. In order to connect the two stages, we take out the profile of a central object 
from the core-collapse simulation and continue the cooling simulation to simulate the shock 
revival and the cease of accretion in the explosion. 

We demonstrated the method to extract the physics information from neutrino emissions 
of the Galactic supernova based on the supernova neutrino database and additional simula- 
tions [376]. We explored the time profile of neutrino detection events at the SK for the set of 
progenitors with different masses and shock revival time. We provided the basic feature of 
the rise and fall of neutrino burst to extract the information of the core bounce of massive 
stars in the early phase within 0.3s. We explored also the long term behavior of neutrino 
burst over 20s up to the phase of fading out. We performed additional long simulations 
of the proto-neutron star cooling in the case of light and massive neutron stars. We found 
that the SK and Hyper-Kamiokande are able to detect the long term evolution over 100s 
for the massive neutron stars in Galactic supernovae. We proposed a method to determine 
the mass of central object by plotting the cumulative event number summed from the last 
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event detection in time-backward manner. These results are used as the basis to prepare for 
the future detection of the next supernova in order to determine the details of supernova 
mechanism and the neutron star. The dependence of the neutrino signals on the equation 
of state (EOS) is now in progress under the development of the EOS table for supernova 
simulations as described below. 


4.2.5. Information of supernova mechanism from supernova neutrino detectors. The accu- 
mulation of predicted signals of supernova neutrinos will enable us to extract the information 
of supernovae from the detection of neutrino bursts in the future [349, 350]. SK will detect 
~10,000 neutrino events from a supernovae in our Galaxy and provide valuable and unprece- 
dentedly detailed information regarding the supernova mechanism. For nearby explosions, 
pre-supernova neutrinos from the silicon burning stage just before explosion will be studied 
by both KamLAND and the newly gadolinium-loaded SK. 

Events from the dynamical phase will reveal hydrodynamical instabilities and rotation, 
which are augmented by the detection at IceCube [317, 351]. Observation of events from 
the cooling phase will clarify the birth of compact object with mass and radius, which may 
constrain the equation of state of dense matter [376]. 

A combined detection of GWs and neutrinos at KAGRA and SK would serve to unveil 
the supernova mechanism by the correlation between them |31, 352]. The planned next 
generation of big detectors such as DUNE, JUNO, and Hyper-Kamiokande [353, 354] with 
their extended supernova neutrino detection ranges can be expected to provide additional 
information (See [355] for further references). 

Regarding multimessenger astronomy, in the case of a galactic supernova an early alert 
is possible since neutrinos are generated earlier in the explosion and therefore arrive before 
the electromagnetic radiation. Such an alert could be generated and disseminated either 
through large, directional, high-confidence detectors like the gadolinium-loaded SK acting 
independently to announce a burst, or a network of detectors which detect lower-confidence 
signals acting in coincidence [355] to reduce the chance of individual false alerts. 


4.2.6. Neutrino bursts from black hole formation. Studying the neutrino burst from the 
black hole formation is an interesting target of the SK among the variety of gravitational 
collapse of massive stars [377]. In the case of stars of more massive than those for ordinary 
supernovae, 40-50Mo. for example, with the intense accretion of matter from outer layers, 
the retreat of shock wave is inevitable and there is no chance to have the explosion. The 
mass of central object increases monotonically and attains the maximum mass, which can 
be supported by the equation of state. The black hole is formed due to the re-collapse of 
proto-neutron star at the critical mass. The neutrino signal in this case has a characteristic 
signature with a short duration, typically about 1s, and increasing average energies and 
luminosities. This is cased by the increasing mass due to the accretion and the associated 
increase of density and temperature. This type of short neutrino burst can be a hallmark of 
black hole formation and can be used to constrain the equation of state of hot and dense 
matter [378, 379]. 

The detection of GWs and neutrino bursts from the black hole formation may provide the 
information on the central object and the equation of state in addition to the case of core- 
collapse supernovae (see sec. 4.1). We revealed the characteristic feature in the frequencies of 
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GWs from the dynamical evolution of proto-neutron star toward the black hole formation. 
We analyzed the time evolution of accreting proto-neutron stars to determine the funda- 
mental and gravity modes of GWs [380] by utilizing the results of the numerical simulations 
of massive stars with a set of equation of state [377]. The density increase toward the black 
hole mass leads to the rise of frequencies of GWs as a function of average density of the 
proto-neutron star. The ratio of the frequencies of the two modes can be characterized by 
the compactness of the proto-neutron star, therefore, the information of mass and radius. 
Moreover, the termination of neutrino signal can provide the timing of the black hole forma- 
tion and the information of maximum density through the combined analysis of GWs and 
neutrino signals. 


4.2.7. Diffuse supernova neutrino background. The numerical simulations of the super- 
nova explosions and the black hole formation from various massive stars are essential for 
the detection of diffuse supernova neutrino background, which is the main target of the SK 
with gadolinium loading. The numerical simulations of the core-collapse supernovae in multi- 
dimensions (see also sec. 4.1) and the black hole formation as discussed above have been 
applied to provide the integrated energy spectra of neutrino emission from the progenitors in 
the wide range of stellar mass. The feature of neutrino burst depends on the compactness of 
massive stars, which in turn determines the accretion rate of matter and affects the neutrino 
emission. In the wide coverage of massive stars with different compactness, the case of black 
hole formation contribute to hard energy spectra due to the energetic bursts. Hence the con- 
tribution of black hole formation for large compactness increases the event rates of diffuse 
supernova neutrino background. The long term observation at SK and Hyper-Kamiokande 
can provide the constraint on the critical compactness and the ratio of black hole formation 
in various stellar collapse [381]. 


4.2.8. Equation of state for supernovae and neutron star mergers. The equation of state is 
one of the important ingredients in core-collapse supernovae as well as neutron star mergers. 
It is also one of remaining uncertainties in nuclear physics to determine the outcome of 
explosions, the birth of compact object and the neutrino emission. In addition to the progress 
of numerical simulations of supernovae, developments of the data table of supernova EOS for 
simulations have been made over the decades. One of the popular sets of supernova EOS is 
the Shen EOS based on the relativistic mean field theory [382-384], which has been applied 
to provide extended versions of data table of supernova EOS [385]. The Shen EOS table is 
publicly available on the web and has set the standard format to provide thermodynamical 
quantities for the usage of numerical simulations. The Shen EOS has been used in the 
numerical simulation for the supernova neutrino database. 

Recently, the Shen EOS has been revised to improve the density dependence of the sym- 
metry energy [386]. This modification is motivated by the recent progress of observational 
data on neutron stars from the GW detection and the X-ray observations. The symmetry 
energy of the original Shen EOS, which was determined by the fitting to the nuclear masses 
and radii, has been claimed to be rather large as compared with other frameworks. Recent 
nuclear experiments also provide constraints on the behavior of symmetry energy and help 
to improve the description of symmetry energy. A new term of density-dependent iso-vector 
interaction in the relativistic mean field theory provides a smaller symmetry energy while 
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keeping good properties of symmetric nuclei and matter. It provides also smaller radii of 
neutron stars within the observational constraints. The updated Shen EOS is now applied 
to numerical simulations of core-collapse supernovae and proto-neutron star cooling. It has 
been recently shown that the influence of the updated Shen EOS at high densities mainly 
appears in the evolution of proto-neutron stars when the matter becomes neutron-rich and 
remains minor in the dynamics around the core bounce. It is interesting to see the influence 
of the full data table of the updated Shen EOS in core-collapse supernovae and neutron star 
mergers. 

In addition to the revisions of the nuclear interaction, the improvement of the supernova 
EOS of hot and dense matter has been made to describe the mixture of nuclei under the 
nuclear statistical equilibrium [387, 388]. There is recent progress of the construction of EOS 
tables by the microscopic nuclear many body theories such as the variational method and 
the Dirac Brueckner Hartree-Fock theory [389, 390] based on the nucleon-nucleon interac- 
tions, being different from the effective nuclear many body theories. In addition, the sets 
of EOS tables with systematic coverage of EOS parameters are under the development and 
will be applied to the study of core-collapse supernovae and neutron star mergers. These 
development will help to study the signal of neutrinos and GWs and to probe the dense 
matter in these astrophysical events. 


4.2.9. Supernovae simulations on GPUs. The dynamics of the supernovae is described by 
hydrodynamic equations for dense matter and the Boltzmann equation for neutrino transport 
under the gravitational effect described by the theory of general relativity. In addition to 
these coupled equations, data of physics processes such as sets of equation of state and 
reaction rates for neutrinos are needed to be implemented in the numerical simulations. 
Since the neutrino distribution is a function in six dimensions (three dimensions in space 
and three dimensions in neutrino momentum space), the required computational resources 
quickly increases with grid resolution of these degrees of freedom. Therefore the numerical 
simulations of supernovae tend to be challenges in large-scale computational techniques and 
algorithms. 

So far most of large-scale simulations have been performed on massively parallel clus- 
ters, such as the K computer. Recently another type of architecture, which makes use of 
arithmetic accelerators such as GPUs, has become popular in high performance computing. 
Although the accelerator device has an advantage in cost performance, the code implemen- 
tation becomes much more involved to achieve desired performance. Furthermore, whether 
accelerators work efficiently depends strongly on the structure of numerical algorithms. 

Presumably for these reasons, application of GPUs to the simulations of core-collapse 
supernovae has been restrictive. As for the simulation code including the neutrino transport, 
the VERTEX code was ported to GPUs by employing CUDA [391]. In the VERTEX code, 
its most time consuming part is calculation of the collision term of the Boltzmann equation 
that exhausts almost half the simulation time. Offloading one dominant reaction process to 
GPUs achieved 1.8 times acceleration of whole simulation time. 

We develop an efficient scheme to exploit accelerator devices such as GPUs in the numerical 
simulations of the core-collapse supernovae [392, 393]. As the first step, we apply the offload- 
ing of simulation bottlenecks to the computations of neutrino-radiation hydrodynamics under 
spherical symmetry [373]. By adopting the implicit scheme for the evolution equation, an 
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iterative linear equation solver for the coefficient matrix is the most time consuming part. To 
offload this part to the GPU devices, we employ OpenACC as a framework of implementation 
as well as make use of cuBLAS library that is available on NVIDIA’s CUDA environment. 
We change the data layout to maximize the efficiency of data transfer between the device 
memory and cores, in accord with so-called coalesced access. With both the OpenACC and 
cuBLAS, significant acceleration is achieved so that the linear equation solver is no longer a 
primary bottleneck [392]. 

The secondary bottlenecks of the original code are the computation of collision term in 
the Boltzmann equation and the inversion of block diagonal part of the coefficient matrix. 
The latter is required in the weighted Jacobi-type preconditioner for the iterative linear 
equation solver [394]. Offloading to GPU devices successfully accelerates these tasks [393]. 
For the inversion of block matrices, we employed a blocked version of the Gauss-Jordan 
algorithms [395] that is suitable to the many-core architecture. With these improvements, 
systematic simulations with better resolution have become feasible. Further optimization 
as well as extension to 2D and 3D simulations are underway. We are also implementing a 
code for PEZY-SC processor, another kind of accelerator device, by applying the technique 
developed for GPUs. 


5. Probes of new physics 

5.1. Test of gravity theories using gravitational wave data 

5.1.1. Varieties of gravitational waveform within general relativity. The GW signals from 
coalescing binaries can be decomposed into inspiral, merger and ringdown waveforms. 

Here we focus on the inspiral phase, in which the orbital frequency and hence the GW 
frequency gradually increases as the binary separation decreases. In the case that the binary 
with aligned spins evolves along a sequence of quasi-circular orbit, the inspiral waveform in 
the frequency domain can be approximately represented as 


A(f) = Af- / 6g) | (13) 


where the overall amplitude A is x M°/®, and 


, (14) 
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in the post-Newtonian (PN) approximation [153]. Here, 7 = u/M is the symmetric mass 
ratio, v = (tM f)!/? is the orbital velocity, and 
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The second term in the right hand side of the above equation becomes a small contribution 
because nyI — 4n < v3/18 = 0.096 - - - , but we should note that this term can be dominant 
for spins antialigned with each other. 

The first term in the square brackets of Eq. (14) is the term obtained by using the New- 
tonian dynamics and the quadrupole energy loss formula due to GW emission. So, when 
we discuss the phase evolution, we count the PN order relative to this leading term. The 
second term is relatively O(v?), which we refer to as 1 PN order. In the following, we sum- 
marize additional effects which should be carefully discriminated from modifications due to 
extension of general relativity (GR). 
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The effect of spins: 
The lowest PN effect of spins on the phase evolution appears at the 1.5 PN order. The 
effect depends only on the spin components projected to the direction parallel to the 
orbital angular momentum at this order, which is encapsulated in the coefficient 5. 
When the spins are not aligned with the direction of the orbital angular momentum, 
spin precession occurs. As a result, the amplitude modulation occurs, which helps to 
solve the degeneracy between spins and other parameters [396]. The GW templates for 
precessing binaries are ready (see, e.g., Ref. [397] and references therein). 

Eccentric binaries: 


19/12 at the late stage of binary evolution, where 


The orbital eccentricity decays like e c a 
a is the orbital semi-major axis [398]. Therefore, binaries are circularized well before 
they merge unless they are borne with a close separation and a large eccentricity. One 
possible scenario to form binaries that maintain a sufficiently large eccentricity at the 
coalescence time is the binary formation by capture in star clusters. Hence, detection 
of eccentricity would be useful to distinguish the formation channels of binaries that 
become GW sources. In the GW waveform, the effect appears most significantly in the 
phase evolution, unless the eccentricity is not extremely large. The correction appears 
in the phase evolution at —19/6 PN order. The minus sign of the PN order reflects the 
fact that the orbital eccentricity decays rapidly. 
Dark matter cloud: 

In the presence of dark matter, the waveform of GWs can be affected [399]. In the 
scenario of the existence of the dark matter mini-spike around an IMBH [400, 401], the 
waveform of the GWs from an EMRI/IMRI system is modified due to the effects of the 
gravitational pull, the dynamical friction, and the accretion of dark matter [402-404]. 
'They have shown a possibility of the measurement of the power-law index of the dark 
matter mini-spike radial profile with space-borne GW detectors such as LISA. (See also 
Ref. [405].) On the other hand, when the mass ratio is around O(10?) ~ O(10*), energy 
dissipation due to the dynamical friction is comparable to the total binding energy 
of dark matter. If one considers a modification of dark matter distribution due to the 
energy dissipation, the differences of the waveform from the vacuum case becomes much 
smaller, even the spike has a large power-law index [406]. 


5.1.2. Model independent test of waveform consistency with general relativity. From the 
inspiral waveform, we can determine the binary parameters. At the same time, we can test if 
there is no deviation from the standard waveform predicted for non-spinning quasi-circular 
binaries in GR. We can test the presence of spin and eccentricity of each binary component 
BH. At the same time we can also test the deviation from the prediction of GR. 

Here we summarize the test of GR performed by LIGO/Virgo collaborations [35]. The 
first test is the one to test if the residual data which is obtained by subtracting the best-fit 
template is consistent with the Gaussian noise. There is no inconsistency observed. 

They also have made a consistency test between the estimations of parameters obtained 
from the inspiral part and from the merger-ring down part. Small perturbation modes of 
a BH have its own frequencies and damping rates determined by its mass and spin. Those 
damped oscillation modes, which are called quasinormal modes (QNMs), are excited at the 
epoch of formation of the remnant BH after coalescence. See, e.g., Ref. [407] as a review of 
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QNMs including the history. Among various QNMs, the dominant mode is the fundamental 
(n = 0) mode with the harmonic (£, m) = (2, 2), and the GW waveform is written as 


h(t) = Ae-*79)/7 cos(2x fa (t — to) — do) , (16) 


where fg is the oscillation frequency, and 7 is the damping time, and tọ and $o are the 
starting time and its phase, respectively. A is the amplitude at the starting time. We can 
find the information of fg and 7 which are related to the mass Mrem and spin arem of the 
remnant BH, in Refs. [408, 409]. 

As for the initial phase, we can maximize the SNR with an analytical method presented 
in Ref. [410, 411] which is based on Ref. [412]. Since it is difficult to determine the starting 
time (e.g., see Ref. [413] for the estimation), the fractional differences 


insp post—insp insp post—insp 
A Mtem M 2( rem — /Mrem ) Aadrem = 2(arem — Qrem ) (17) 
y 7 ins ost—ins ? A D ins ost—ins ? 
Mrem Mren + ME p Qrem Gian 3 afem p 


between the quantities for the inspiral and the post-inspiral phases have been used as an 
inspiral-merger-ringdown consistency test for BBHs in Ref. [35]. In the GR case, we will 
have 


AMrem =f) AO eni 


, - 
Miem Qrem 


= 0, (18) 


within the measurement error range. The analysis result for GW events in the LIGO-Virgo 
Catalog GWTC-1 is shown in Fig. 2 and Table III of Ref. [35]. Here, we should also note that 
the SNR only for the ringdown phase is not so high in the current observations. Therefore, 
the merger-ringdown phase with a combination of phenomenological and analytical BH per- 
turbation theory parameters [414], is treated in the parameterized test of GWs. This analysis 
result for GW events in the LIGO-Virgo Catalog GWTC-1 is also shown in Fig. 3 of Ref. [35]. 
On the other hand, when we estimate the mass and spin of the remnant BH (e.g., the right 
panel of Fig. 4 in Ref. [5]) without testing gravity, we use the inspiral-merger-ringdown 
waveform from numerical relativity [415-417]. 

Using only the ringdown phase, one of the best way to test gravity is “black hole spec- 
troscopy” by observing multiple QNMs [418] (see also Refs. [419, 420] for testing the no-hair 
theorem with black hole ringdowns, and another simple method is proposed by Ref. [421]). 
For example, we may consider Eq. (17) with substitution, 


insp _ L=2,m=2 post—insp _ £=3,m=3 
Mizse — M !, MPS = Men), 


rem rem ? ( 
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if the (€ = 3, m = 3) QNM is the next dominant mode. In order to do so, it is necessary to 
extract fg and 7 of the weak ringdown signal from the noisy data accurately for each QNM. 
In Ref. [422]. As the first step, mock data of GWs which include some deviation from the 
GR prediction was prepared, and it was applied to the following five methods to extract 
the dominant QNM: (1) plain matched filtering with ringdown part method, (2) matched 
filtering with both merger and ringdown parts method, (3) Hilbert-Huang transformation 
method, (4) autoregressive modeling method, and (5) neural network method. It was found 
that the determination of fg is much easier than that of 7 although the accuracy depends on 
the analysis methods. Interestingly, it turned out that the standard matched filtering does 
not give the optimal inference in this problem. 
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As a unified framework to test possible modifications of gravity, one can use parametrized 
post-Einstein (PPE) approach [423]. In this approach, we consider the modification of the 
functions that appear in the GW waveform in GR (13) without specifying the model as 


(1+ > o, )Agn(f), Y) > var(f) +X Biv”, (20) 


where q; and 8; specify the amplitude of modification, while a; and b; specify the PN order of 
modification. The point is that the leading order correction tends to be given by a particular 
PN order. In many cases the leading order effects appear in the phase evolution. Of course, 
the test targeting at a particular modification can be more sensitive. 

LIGO/Virgo collaborations did not give the constraints on the PPE parameters directly. 
Instead, they examined the constraint on the modification to the model parameters contained 
in the waveform model IMRPHENOMPv2, following the method proposed in Ref.[424] and 
implemented by Ref.[425]. For the inspiral part, they test the modification to the existing 
PN coefficients in the phase evolution up to 3.5 PN order in GR as well as the —1PN and 
0.5 PN coefficients which are not present in GR templates. Also, the modification of the 
parameters that characterize the merger and ring-down phases is also tested. 

For the modification of GW propagation, what they did is similar to PPE, but the con- 
straints are derived for the amplitude of modification per propagation distance by stacking 
the data. In all cases mentioned above, significant deviation from GR prediction has not 
been obtained. 


5.1.3. Model independent test of extra polarizations. 'The GW interferometers basically 
detect the tidal deformation of the physical distances among nearby test masses, which can 
be described by the geodesic deviation equation for slowly moving objects 


i vi —R' jor! » (21) 


Namely, what we can observe is the rank 2 tidal tensor Ej; = jojo. The absolute distance 
cannot be measured by the current GW interferometer, and they measure the difference 
of the distance changes in different directions. Hence, the trace part of E;; is irrelevant. In 
generic dynamics of geometry, therefore, there are five components which can be decomposed 
into three dimensional 1 scalar, 2 transverse vector, 2 transverse traceless tensor components. 

It is not easy to think of gravity models in which scalar or vector components of the tidal 
waves are emitted without modifying the GW waveform significantly. For example, in the 
case of scalar-tensor theories, the effect of scalar wave emission, which also contributes to 
the excitation of the scalar component of the tidal tensor, appears most significantly in the 
orbital evolution of the binary due to the radiation reaction that it causes. If we consider 
such models that can excite larger magnitude of tidal tensor in the scalar or vector mode 
without losing much energy, the coupling of matter field to those modes tends to be too 
large to hide the influence on the tests of gravity in solar system. 

For the reason mentioned above, sensible models would be a conversion of the excitation 
energy from the tensor modes, i.e., standard GW modes, to the other exotic modes during 
the propagation. In that case, the waveform can be identical to the original GW waveform 
although the phase velocity of the scalar and vector waves can be different from the speed of 
light. In that case, it is difficult to put some constraint to the conversion of the energy during 
the propagation. So far, only comparison of likeliness between pure tensor modes and pure 
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scalar mode or between pure tensor modes and pure vector modes has been done, finding 
that pure tensor modes are largely preferred [35, 426, 427]. 

In a general framework, a GW signal can be a superposition of scalar and/or vector modes 
in addition to the ordinary tensor modes: 


St.) = S Nr ATO. (22) 


Here FP is the antenna pattern function of the a-th detector for the polarization P, and 
hp is the GW waveform for the polarization P, where P= +, x, V, W, S and L [428-430]. 
V and W modes are vector modes and S and L modes are scalar modes, which satisfy 
F? = —FL [428]. The antenna pattern functions are dependent on a GW source location 0 
and 9, which are the latitude and longitude, respectively. The separability of the polarization 
modes of GWs from a point source has been investigated by solving an inverse problem for the 
first time in [431]. It has been shown that in a nontensorial polarization search, the necessary 
number of detectors is at least the same as the number of polarization modes to be searched, 
e.g., at least three detectors are necessary to search for a superposition of two tensor and 
one scalar modes. The mode separability in the case of compact binary coalescences is subtle 
because their GW waveforms are well modeled with multiple source parameters. It is shown 
in [432] that even with correlations and degeneracies among the parameters, a mixture of 
the polarization modes are separable with the same number of detectors as the number of 
the modes. Then the scalar or vector modes with similar amplitudes to the tensor modes can 
be detected. For a binary neutron stars (BNS) observed with the third-generation detectors 
such as Einstein telescope[433] or cosmic explorer[434], the mode separation with smaller 
number of detectors is possible due to the long duration of a signal (Z 1 hour) and the 
time-evolving detector antenna patterns [435]. 

With the help of an electromagnetic counterpart, more polarization modes than mentioned 
above can be searched once the sky location is fixed. We will discuss this possibility in the 


next. 


Direct polarization test with the help of an electromagnetic counterpart. We consider 
to analyse the data from a GW source, S,(t,Q) + nalt), where nq is the detector noise. 
Equation (22) means that tests of all polarization modes with a network of interferometers 
require at least five detectors. However, Hagihara et al. pointed out that there exist particular 
sky positions that allow for a vector mode test, because scalar modes, even if they exist, 
can be perfectly eliminated from a certain combination of the strain outputs at the four 
detectors only for these sky directions [436] (See Fig. 35). They also found that a vector 
mode test can be done by using four detectors [437]. They put also a direct upper bound on 
a vector component from GW170817 event, since the scalar modes can be largely suppressed 
for this event [437, 438]. 

Very recently, they proposed a general formulation for directly testing extra GW polar- 
ization and a certain condition that allows for a scalar test. First, we focus on the network 
composed of the four ground-based interferometers. For a given source, we know its sky posi- 
tion, as is the case of GW events with an electromagnetic counterpart such as GW170817. 
'Then, one knows exactly how to shift the arrival time of the GW from detector to detector. 
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Fig. 35 Seventy sky positions in the earth-centered coordinates. At these points for a GW 
source, the spin-1 test can be done in principle, because the spin-0 and spin-2 modes can be 
perfectly eliminated. The numerical calculations were done in Reference [436]. 


For example, the projection operator is defined to eliminate W, 4- and x modes. It 
is denoted as II?W**, For this example, hw, h} and hy in the strain outputs {S,} are 


eliminated as 
MWH S, = (ettet SEI ET FF) (hg — hy) + (E EY FI RS FT) hy HIY na, (23) 


If the coefficient of hy in Eq. (23) vanishes in a certain sky region, there remains only the 
spin-0 part in the null steam. Therefore, the spin-0 polarization test is possible, if à GW 
source is found in this sky region. The vanishing coefficient condition is 


FY FW Get ls 
FY FY o F 
BY FY gu FZ 
Fy Fy Fi Fg 


=0. (24) 


Note that this is invariant for choosing a reference axis for the polarization angle. D4 = 0 (or 
sufficiently small D4 practically) is a condition for directly testing scalar modes separately 
from the other modes only by four detectors [438]. 

A possible straightforward procedure of the GW data analysis along this direction is as 
follows. At first, we determine the sky location of the GW/EM source from multimessenger 
astronomy, especially by optical and VLBI observations. Secondly, by using the sky location, 
the arrival time shift for each GW detector is taken into account. Next, the (time-shifted) 
strain output at each detector is substituted into Sa in the left-hand side of Eq. (23). If the 
left-hand side of Eq. (23) is above the noise level (the last term of the right-hand side), then, 
the existence of extra GW polarizations could be suggested. If it is comparable to (or lower 
than) the noise level, a certain direct upper bound can be placed on extra GW polarizations 
(through the first and second terms the right-hand side). In this procedure, explicit templates 
of GW waveforms are not needed. In this sense, this test of extra GW polarizations is robust. 
Further study along this course is left for future. For instance, sophisticated algorithms and 
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Fig. 36 Contour plot of D, in the earth-centered coordinates. The numerical calculations 
were done in Reference [438]. 


pipelines would be important for real data analysis, because the nature of detector noises is 
quite complicated. 

With four detectors, consistency tests with GR can be done by examining if the wave in the 
fourth detector, say KAGRA, is consistent with the three-detector network. Equation (22) 
for strain outputs (a = 1,--- ,5) from the five detectors including LIGO-India can be always 
solved for hs — hr, hy, hw, hy and hx. In principle, adding the fifth detector will thus allow 
for a direct search of extra GW polarizations for any sky region. This would be interesting 
in probing new gravitational physics beyond GR. 


5.1.4. Gravitational waveform corresponding to various modifications of gravity. The 
sources of GWs are composed of the objects of extremely strong gravity in the sense that 
the deviation from the Minkowski spacetime is significant. The spacetime curvature radius 
is comparable to the size of the system in the case of BH and/or NS binaries. 

In the case of BBHs, the emission of electromagnetic counterpart is not expected much. 
In fact, no detection of electromagnetic counterpart has been reported so far. In the case of 
binaries that include a NS, we have a chance to detect the electromagnetic counterpart. In 
fact, various follow-up observations succeeded in finding the counterpart for GW170817, the 
first coalescing BNS event detected by GWs [9] (see Sec. 3). However, even in that case the 
emission of the electromagnetic counterpart in the region of truly strong gravity cannot be 
observed directly because of the large opacity of the high density matter. 

In contrast, GWs are emitted from the bulk motion of the objects where the strong grav- 
ity is really at work. Hence, GWs are thought to be an indispensable probe of strong 
gravity. Here, we summarize various modifications on the GW waveform by considering 
representative extended theories of gravity. 


Scalar tensor gravity. As the most typical modification of gravity, one can consider scalar- 
tensor theories. A representative theory will be Brans-Dicke theory [439], defined by the 
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action 
1 


Se ram dryg (¢R- wend $ ap”) - Y | sto» (25) 


where wgp is the so-called Brans-Dicke parameter. When we send wgp to oo, Einstein 
gravity is recovered. m,(¢@) is the effective mass of a star labeled by a, and it depends on the 
local value of the scalar field around the star. The effective gravitational coupling depends 
on wep as Gef = (4+ 2wpp)/($(3 + 2wpp)), and the dimensionless scalar charge is defined 
by s, := Olnm,(¢)/Oln ¢, which can be interpreted as the dependence of the gravitational 
binding energy of the star on Gog. 

The modification to the GW waveform appears at —1PN order like (f) =---+ 
(3/128) (rM f)>/8 [o 72/3 +1+---] in Eq. (14), where the coefficient a is given by a= 
—5(s1 — s3)?/(64wpp) [440]. This is caused by the dipole radiation of the scalar field. If the 
scalar charge of two stars are proportional to their mass, i.e., if the dimensionless scalar 
charges are identical, the dipole radiation does not occur. Hence, the dipole radiation is not 
expected from equal-mass BNSs. 

In this model, numerical approach is also possible, but the effect is larger for lower fre- 
quencies. Therefore, the template based on the PN approach mentioned above would be 
sufficient, unless positive detection of the presence of scalar charge is reported. 

One scenario to give a scalar charge to compact stars is spontaneous scalarization in the 
context of generalized Brans-Dicke model, in which wgp is promoted to a function of $ [441]. 
Rather intuitive understanding of this process can be obtained by introducing the canonically 
normalized field y = f \/2wpp(¢)/¢d¢. In terms of this new field, the gravitational action 
would be rewritten as 


S= is; | tv (or - 0.00") (26) 


From this expression, one can find that the y-dependent part of the energy of a static body 
would be given approximately by 


f dx (pap — 8roly)p) ~ R? (p/ R? — 8rpó(o)) , (27) 


where, after the curvature is replaced with the energy density p using the Einstein equation 
as an approximation, the deviation from the Minkowski metric is neglected. In the last line, 
the size of the system is set to R and q is assumed to vanish outside the compact object. If 
the function ¢(y) is convex upward at the origin and the density becomes high enough or 
the size of the compact object is large enough, non-trivial configuration of the scalar field 
develops inside the star, which leads to giving a scalar charge to the star. The point is that 
this happens only in the situation where pR? zz M/R is sufficiently large. 

Within the category of scalar-tensor theory, more generalization is possible. The Horndeski 
theory [442] is the most general scalar-tensor theory having second-order equations of motion. 
This theory can, however, be generalized further in a healthy way to the so-called degenerate 
higher-order scalar-tensor (DHOST) theories [443-445], which inevitably have higher-order 
equations of motion, but still propagate the same number of degrees of freedom as the 
Horndeski theory does (i.e., one scalar and two tensor modes). The Horndeski and DHOST 
theories offer us a powerful unifying framework to study the physics of dark energy and 
modified gravity (see, e.g., [446] for a review). 
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Measuring the speed of GWs can constrain the Horndeski and DHOST theories. The 
Lagrangian for the DHOST theories satisfying caw = 1 is given by [37-40] 


L = Go(¢, X) — G3(¢, X)Llo + f(b, X)R + Az(¢, X Opo pup” 
+ Aappu” by + As(b4 dur ")?, (28) 


where X :— —o"¢,,/2 and 


1 A 
-35 (2Asf - 3f& — 24s fx X + ASX?) , As = — ^7 (fx + AX). (29) 
We thus have 4 free functions of ¢ and X. For A3 — 0 and f = f(9), this reduces to a 
subclass of the Horndeski theory. Solar System and astrophysical experiments/observations 


Ag = 


are potentially able to put further bounds on the functions in the above Lagrangian. 

For Solar System and astrophysical tests of the Horndeski and DHOST theories, it 
is important to understand the Vainshtein screening mechanism [458], which suppresses 
the scalar-mediated force and typically is implemented in scalar-tensor theories whose 
Lagrangian depends on the second derivatives of the scalar field. In the Horndeski sub- 
class with A3 = 0 and f = f($), the Vainshtein mechanism is known to work perfectly for 
a static spherically symmetric body, ® = V = —Gy M/r, inside a certain radius (which is 
sufficiently larger than the size of stars and the Solar System), where ® and W are two metric 
potentials defined by 0go9 = —26 and ógi; = —2WÓjj, Gn :— 1/167 f is the effective Newton 
constant, and M is the mass contained within r. Therefore, it would be difficult to place 
Solar System/astrophysical constraints on this subclass. 

In generic DHOST theories, the behavior of weak gravitational fields is more interesting. 
It was found in [41-43] that 


M 5Y: M' 
/ 7, / n 
v= 6r (ear), W =Gy(- tT ena), (30) 


where Gy and Y; are written in terms of f and A3, and the prime denotes differentiation 
with respect to r. Since M — constant outside the source body, Eq. (30) shows that gravity 
is modified only inside astrophysical bodies in generic DHOST theories. (In other words, 
Vainshtein screening operates nicely in the vacuum exterior region.) This partial breaking 
of Vainshtein screening was first discovered in [447]. A number of astrophysical tests of 
Vainshtein-breaking theories have been proposed so far. Going beyond the weak gravity 
regime, it was shown in [448] that the structure of relativistic stars is quite sensitive to 
the Vainshtein-breaking effect, which would potentially give tight constraints on DHOST 
theories in the strong gravity regime. 

It was pointed out in [449, 450] that gravitons can decay into ¢ in DHOST theories. To 
avoid this, it is required that A3 — 0. Otherwise, GWs would not be observed. Upon imposing 
this constraint we have A4 = 3f%/2f and As = 0. This yields a special subclass of DHOST 
theories characterized by 3 free functions, in which the behavior of weak gravitational fields 
is completely different from (30), as shown in [451, 452]. (In fact, Eq. (30) is derived assuming 
As #0.) The main result of [451, 452] is summarized as follows: (i) a kind of fine-tuning is 
required even in the vacuum exterior region so that Solar System tests are evaded; (ii) in the 
interior region the metric potentials ® and V obey the standard inverse power law, but the 
two do not coincide; (iii) the interior and exterior values of the effective Newton constant 
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are different, 
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All of these features are in sharp contrast with the case of generic DHOST theories. 


Gin = (31) 


A stringent bound on this particular class of DHOST theories with A3 = 0 comes from 
the fact that the measured value of the Newton constant (which is supposed to be Gext) is 
different from the gravitational coupling of GWs by a factor of 1 — X fx /f. This can be seen 
from the quadratic Lagrangian for GWs, 


_, 
~ 1670 f° 


Oy = eee (MMO be) + PP Tw Gow (32) 
Assuming that the scalar radiation does not take part in the energy loss, one may naively 
replace the Newton constant in the standard quadrupole formula with Gow. Then, from 
the orbital decay of the Hulse-Taylor pulsar we find that |X fx /f| € 107%. It would be 
interesting to explore other astrophysical constraints on this special class of Vainshtein- 
breaking modified gravity. 


Quadratic gravity. As a low energy effective theory of the dynamics of the metric, the 
term in Lagrangian at the lowest order of mass dimension is the cosmological constant and 
the next leading order is the Einstein-Hilbert action. In this sense, possible natural extension 
of gravity is to consider to add terms quadratic in curvature to the Lagrangian. 

There are three parity preserving terms, R?, Ruy RY and Ryypo R”. R? term is classified 
as an alternative representation of scalar-tensor theories [453]. One combination, Reg = 
R? — 4R R"" + Ruvpo R}”®? , is know as the Gauss-Bonnet term, which does not contribute 
to the equations of motion as it can be written as a total divergence. Hence, truly remaining 
new possibility is, e.g., adding Ruy R"". However, this theory has a ghost degree of freedom, 
i.e., it contains infinitely many negative norm/energy states, which leads to an absence of 
stable vacua. This is caused by the presence of higher derivative terms in R,,, R"". Hence, 
we need to treat such a model simply as a low energy effective theory. Then, the higher 
derivative terms can be replaced by using the lower order equations of motion, which is 
now the Einstein equations, up to an appropriate field redefinition [454]. Once we apply this 
procedure to replace Ry R"" with T,,,T"", the model is reduced to just Einstein gravity 
with matter fields described by a little complicated Lagrangian. 

As a parity violating term quadratic in the curvature, we can think of the Chern-Simons 
term, RR :— eP V ROX pv RY” ,5/2, where &""?? is the Levi-Civitá tensor. However, this term 
is also written as a total divergence, and hence it does not contribute to the equations of 
motion. 

One interesting extension is to introduce a scalar field coupled to these higher order cur- 
vature terms. The two cases with a scalar field coupling to the Gauss-Bonnet term and the 
Chern-Simons term have been extensively studied. The Lagrangian is given by 

1 rs 


E l ecu 2 
~ 167Gn R+ 7 (Rep or RR) - (06) —V(9)| . (33) 


where £ is the length scale that characterizes the coupling of the axion field to the Chern- 
Simons term. The deviation from GR in such theories can be easily hidden in the weak 
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gravity regime, An interesting phenomena in this setup is that BHs can have a scalar hairs. 
The constraint coming from the GW observations is discussed in Sec. 5.1.5. 


Massive gravity. If we simply add a mass to the graviton, the equations for the GW prop- 


agation in flat spacetime Ohy = 0 would be expected to be modified to (O — m?)huy = 0. 
However, this simple-minded extension is not theoretically very attractive. One consistent 
way that preserves the Lorentz symmetry on Minkowski background is the dRGT ghost-free 
massive gravity [455]. However, it turned out to be difficult to construct even the homoge- 
neous isotropic universe model in this framework, and hence various generalizations were 
proposed. In any case, the dRGT ghost-free massive gravity model requires auxiliary metric 
or tetrad that do not transform as dynamical fields. As a result, the general covariance is 
not maintained. 

In this class of models of massive gravity, linear perturbation predict a large excitation of 
a scalar graviton around a star, even if we send the graviton mass to a small value, which is 
known as the van Dam-Veltman-Zakharov discontinuity [456, 457]. However, in this limit, 
the linear analysis cannot be trusted any further and the non-linear effects become large, 
strongly suppressing the excitation of the scalar graviton. This is the very same mechanism 
as the Vainshtein mechanism introduced above [458]. As a result, the gravity around massive 
objects is thought to behave very similarly to GR. 


Bigravity. A natural way to recover the general covariance in the theory with massive 
graviton is to promote the auxiliary metric field in the dRGT ghost-free massive gravity to 
a dynamical one. Such a model is know as dRGT bi-gravity theory, which is also shown to 
be ghost free by Hassan and Rosen [459]. However, the general covariance guarantees the 
presence of a massless graviton. In this model the field that becomes massive is the second 
graviton, which corresponds to the difference of two metrices, roughly speaking. 

The generation and propagation of GWs from binaries were studied in Ref. [460], in which 
the possibility of graviton oscillation between two graviton modes: one is massless while the 
other is massive, was pointed out. In the ordinary situation, only the metric directly coupled 
to the matter field is curved, while the other remains to be almost flat. In a particular choice 
of the model parameters, however, the massive graviton becomes very heavy around a star, 
and hence the difference between the two metrices is suppressed. In this case, the graviton 
oscillation becomes possible, without violating the test of gravity, say, tests in the solar 
system. 


Einstein-aether theory. The Einstein-aether theory, which for hypersurface-orthogonal 
configurations of the aether vector can be considered as the low energy limit of a candidate 
theory of quantum gravity called non-projectable Horava-Lifshitz gravity and which also 
serves as a testing ground of GR. We thus review the theoretical and observational constraints 
on the Einstein-aether theory The basic variables of the gravity sector are [461] the four- 
dimensional metric gv, the aether vector u” and a Lagrange multiplier A, where u,v = 
0,---,3 and we adopt the signatures (—,+,+,+) for the metric. The total action is [462] 
S = Se + Sm, where 


Se = age | VT te [Riow) + £e (gu) ]; Sm = f V58 dafen ton ]. G9 
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Here v denotes matter fields, R and g are the Ricci scalar and the determinant of gav, and 


Le = —M „y (Dau) (Dgu^) +A (gagu*u? +1) , (35) 

where D,, denotes the covariant derivative compatible with guv, and 
Me wo = cig”? guv + 02506) + czo oh — c4u* uPg,,, . (36) 
It is convenient to introduce linear combinations of c;'s (i = 1,...,4) as cij = cj +c; and 


Cijk = Ci + cj + cy. The Newton constant G y is related to Ge as [463] Gy = G4/(1 — 5C14). 

The Minkowski spacetime with u, = on is a solution of the Einstein-aether theory. It is 
then straightforward to analyze linear perturbations around the Minkowski background and 
investigate properties of spin-0, -1 and -2 modes. The coefficients of the time kinetic term 
of each excitation must be positive, 


1 — c13) (2+ c13 + 3C2 
qs vr > 0; eas — ) qv —C4, Qr=1-—ci3. (37) 


In addition to the ghost-free condition, we must also require the absence of gradient insta- 
bilities by demanding that the squared speeds of propagation be non-negative. Moreover, 
the almost simultaneous detection of GW170817 [9] and GRB 1708174 [36] sets a strin- 
gent constraint on the speed of the spin-2 mode as —3 x 107 < ep — 1 < 7 x 10716, where 
c? = 1/(1 — c3). This implies 

lag] < 10715. (38) 


As for spin-0 and -1 modes, the squared speeds must be greater than 1 — O(1) x 10715, in 

order to avoid vacuum gravi-Čerenkov radiation [464]. We thus impose 

2c1 — 13 (26i = c13) 
2c14(1 T C13) 


c123(2 — c14) 
ci4(1 — c13)(2 + c13 + 3e2) ’ 


ye , (89) 


; 2r. 
Cov Z+; Cg = 


in addition to (38). 
Applying the theory to cosmology, one can find that the gravitational constant appearing 
in the effective Friedman equation is given by [463], 


1 
Ge p= 314 


sca — = 40 
aa $ (cig + 3e2) 1+ GE + 3c2) o 
One then needs to impose the nucleosynthesis bound as 
Geos 1 
—— —1' <S. 41 
35 1$ (41) 


Among the 10 parameterized post-Newtonian (PPN) parameters [465], the only two 
parameters that deviate from GR in this context are 


8(c$ 1 2c3 — ca)(2c, +3 
ese (G T oa) opes (c1 + 2c3 — c4) (2c1 + 3e + c3 + c4) (42) 
26; — cÍ + c$ 2 c123(2 — c14) 
In the weak-field regime, Solar System observations set the constraints [465] 
lo1| < 107%,  |os| € 1077. (43) 


In the strong-field regime, the isolated millisecond pulsars PSR, B1937 + 21 [466] and PSR 
J17441134 [467] set the constraint [468, 469] as 


lài| < 107°, |â2| < 10? (95% confidence) , (44) 
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Fig. 37 This figure, adopted from [473], shows the constraints in the (co, c14)-plane 
summarized in (47) with two different scales. 


where (@1, @2) denotes the strong-field generalization of (o, 02) [470] given by [471], 


om (8 + a1)Ce . ; àj—04 (c14 — 2)(o4 — 202)0 
2 2(c14 = 2c13) 


ài — 014 , 2-024 
Here, c is the sensitivity. Unfortunately, the sensitivities o, of a neutron star, which depend 


T (45) 


on c;’s and the equation of state of nuclear matter [471], are not known for the new ranges 
of the parameters. We thus leave the analysis of these two constraints to a future work. 
Putting all theoretical and observational constraints together, we have (38), as well as 


0 < c14 < 2.5 x 1075, C14 SZ C1 (46) 
in the (c1, c14)-plane and 


C14 (C14 + 2c2C14 — c2) 


0< c14 Sa < 0.095, —1077 < 
C2 (2 = C14) 


< 1077 (47) 


in the (c5, c14)-plane. In Fig. 37, we show the constraints (47). All the constraints except 
(38) are divided into two classes, those in the (c1, c14)-plane and those in the (cs, c14)-plane. 
For further discussions about the constraints such as the comparison with those on the 
hypersurface-orthogonal theory [472], see [473]. 


5.1.5. Dedicated test for particular models. To optimally constrain particular models, it 
would be better to employ dedicated test. One can imagine various extensions of the simple 
GW waveform adopted in the PPE framework. 


Massive dipolar radiation. One simple example is to add a mass to the field that extracts 
energy from the system through the dipole radiation. Since the radiation is emitted only in 
the frequency range above the mass scale, the modifications in the orbital evolution and hence 
in the GW phase evolution show up only at high frequencies [474, 475]. Therefore, even if 
a NS has a scalar hair, such scalar-tensor theories can evade the constraint from pulsar 
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observations due to the mass of the scalar field. A simple model of scalar-tensor theories 
with a massive scalar field has been constrained by several observations: (1) m € 107?"eV 
for stability on cosmological scales [476], (2) m < 107 !6eV from the observations of a pulsar- 
white dwarf binary [477], and (3) 107 PeV X m € 10-!!eV [478, 479], which relies on the 
measurements of high spins of stellar mass BHs [480]. 

As mentioned earlier, Einstein-dilaton-Gauss-Bonnet (EdGB) theory, is an interesting can- 
didate of modification of gravity obtained by adding a quadratic term in curvature. We need 
to introduce a scalar field, which we call dilaton here, to make the Gauss-Bonnet term rele- 
vant. In this model a BH can possess a scalar charge and hence dipole radiation is expected. 
'The constraint from observed GW data on the model in which the scalar field is massive was 
examined in Ref. [481]. In the massless limit, on the other hand, the current constraint on 
EdGB theory is /agaap Á 1.9 km, which is obtained by using low-mass X-ray binaries [482]. 

Activation of dipole radiation due to massive fields and modifications to GW waveform has 
been analyzed for the GW events in the catalog GW'TC-1 [483] to constrain the magnitude 
of the modification nd [481]. Since the ground based detectors, LIGO, Virgo, and KAGRA, 
are sensitive at frequencies around 10-1000 Hz, GW signals for those are useful to test dipole 
radiation of the massive field in the range 10-!^eV X m € 10-7 ?eV. For m € 10-4^eV the 
modification of GWs effectively reduces to the massless limit. 

'Taking into account the dominant correction to the energy flux due to the dipole radiation, 
which appears at —1PN order, the GW waveform in the frequency domain, h( f), can be 
expressed as [484, 485] 


h(f) = han + &A) &?*. (48) 
The corrections to the amplitude and the phase are, respectively, expressed as [481] 
1 i? —1 3/2 a 7 
ðA = A u ? O(a? — 1), (49) 
5 j AS 3 
zc cer 3) @(@? —1 
ow 18 (a) F(@) O(@ ), (50) 
with 
AA cn ES op a Gy ^12 3/2 


where w = m f/m, m is the mass of the field, O is the Heaviside step function, A is a parameter 
denoting the amplitude of the dipole radiation relative to the quadrupole one. The modified 
waveform, Eq. (48) with Eqs. (49) and (50), can model the modification due to the vector 
field as well as the scalar one. In this parametrization, there are two free parameters, i.e., 
the mass of the field m and the relative amplitude of dipole radiation A. Note that A is a 
function of the system parameters, such as the masses and spins, and the parameters of the 
modified gravity theory that one considers. 

Using the above waveform as a template, 9096 confidence level (CL) constraints on A 
for each LVC GW event has been derived [481]. Since the modified waveform is valid only 
for the inspiral phase, the calculations were terminated when the mass scale reaches an 
appropriately chosen threshold frequency for each event. We should recall that in general 
a strong constraint on A does not always imply a strong constraint on the modification to 
theory. This is because A is proportional to the squared difference of the scalar charges of 
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constituents of the binary as we will see in Eq. (52). Therefore, GW events with vanishing of 
this difference are not sensitive to the modification, even for a quite large coupling constant. 
The constraints on A can be translated to a constraint on the coupling strength, using the 
relation [486—488] 
4g 57 Oba (m3sB4C — uu - 
6 M4 Mq18/ 5 2 


where 


Jas 2V1- X — q.15) (53) 
dg com : 

describes the spin-dependence of the BH's scalar charges. After combining five BBH 
events possessing relatively high SNR, i.e., GW150914, GW170814, GW170608, GW170104, 
GW151226. J/ogaap is constrained as J/agaap S; 2.47 km for all mass below 6 x 10- MeV 
for the first time with 9076 CL including \/agagp S 1.85 km in the massless limit. This con- 
straint in the massless limit is much stronger than the results of in Ref. [489, 490], while it 
is accidentally almost the same with that by low-mass X-ray binaries [482]. 

Since the number of events with a lighter chirp mass, such as BNSs, must increase in 
the near future, it is interesting to stack those events by assuming various theories in 
which charged NSs are allowed consistently. Furthermore, the NS-BH binary is also interest- 
ing [491], because, in addition to their light chirp mass, A can be large in a theory such that 
only either NS or BH has a charge, as in the case of EdGB theory. Moreover, future multiband 
observations will improve the current upper bounds of various modified theories [492]. 


Propagation of GWs in parity violating gravity. The nearly simultaneous detection of 
GWs and the gamma-ray burst from the merger of neutron stars, GW170817/GRB170817 
A [9, 36], offered us an unprecedented opportunity to measure the speed of GWs, caw, at the 
level of one part in 10!°. This puts a tight constraint on modified gravity as an alternative 
to dark energy. We now discuss the impact of this constraint on parity-violating sector of 
gravity [493] 

Let us consider parity-violating gravity whose action is of the form 


| 1l 
^ 167G 


f PE RAD (54) 


where £py is the parity-violating Lagrangian and Ly is the Lagrangian for a scalar field 
@ which may be coupled nonminimally to gravity. The most frequently studied example is 
Chern-Simons (CS) gravity, for which Lpy is given by 


Lev = Los := f (0) RR, (55) 


Recently, it was found that one can construct other parity-violating Lagrangians for the 
gravity sector such as [494] 


4 
Ley = X _arlh, 649") LA, (56) 


I=1 
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with Ly := eve Ro iru X0 ^, Lo := eo Rape Tu buds L3 := purae Ragoo R? OP Pu, 
L4 := ¢)¢*P, and oy :— Vuo. If 4a, + 2a? + a3 + 8a4 = 0 is satisfied, the dangerous Ostro- 
gradsky ghosts are removed in the unitary gauge.” 

We consider GWs hj;(t, T) propagating on a cosmological background. The above examples 
yield the quadratic Lagrangian of the form 


2) _ 1 [a(é) 
a 


M t 
£9) 1 x Ah; Ahr + BC ) 


c3 € P hüdjhii ; (57) 


where €^ is the antisymmetric symbol, A is some energy scale, and a and B are dimensionless 
functions of time determined essentially by the cosmological evolution of the scalar field. In 
the most familiar case of CS gravity, we identically have a = 8. However, we may have 
a # B in general. The Lagrangian (57) can thus be used to describe the GW propagation 
in parity-violating gravity in a unifying way [493] (see also [496] for an effective-field-theory 
approach to arrive at the same Lagrangian). 

In the presence of the parity-violating interactions (57), the + and x polarization modes 
are coupled in the equations of motion, but one can obtain two decoupled equations by 
employing the circular polarization basis. In general parity-violating gravity, the equations 
of motion for GWs in Fourier space can be written as 


(ho)" + (2+ vAYH(h2 + (cáw) K^? hà = 0, (58) 
where the prime denotes differentiation with respect to the conformal time, H :— a’/a, and 


uc Aala — d'H!) k (Ay)? = 1— AABk/a 
4^ ^1-Aa4ak/aÀ ak’? “CW ~ 1—Xaak/aA ’ 


Here, A labels left and right circular polarization modes and A4 = +1. It can be seen that 


(59) 


CS gravity corresponds to the special case with (ras) = 1. Therefore, from the propagation 
speed of gravity no constraint can be imposed on CS gravity. However, parity violation in the 
gravity sector in general leads to (ele # 1 and thus is tightly constrained by GW170817. 
For k/a ~ 100 Hz, one has 


Aa = 8| € 10- km. (60) 


This is the new limit on parity violation in the gravity sector derived in [493]. 

By parametrizing the modification of GW waveforms during propagation in this theory, 
the data in the 01/02 catalog has been reanalyzed by a grid search with the matched 
filtering method [497]. The results imply that the above constraint improves by 7 digits 
since the matched filtering is sensitive to the phase modifications, while the constraints on 
the amplitude modification are similar to the previous result by Yagi and Yang [498]. 

For a face-on binary, the effects of parity violation of gravity cannot be constrained because 
only one of the circular polarization modes can be observed. Therefore, obtaining a more 
stringent constraint or find the signature of parity-violating gravity requires the detection 
of GW signals from a nearly edge-on binary. 


? The ghost degrees of freedom might not be problematic if the theory is treated as a low-energy 
truncation of a fundamental theory. Although these ghost modes do cause instabilities if the theory is 
regarded as a complete theory, it is argued in [495] that the would-be ghost mode that is not manifest 
in the unitary gauge may be an instantaneous mode which does not in fact propagate and hence is 
not dangerous. 
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Echoes from merger remnant black holes. Exotic compact objects are proposed as alter- 
native models to BHs. They are assumed to be compact enough to possess the light ring, 
but without the event horizon. In such models, the feature of GWs after the ringdown 
phase of compact binary coalescences is considered to be different from the BH case. The 
infalling waves generated during merger-ringdown phase might be reflected at the surface of 
the object, which would be just absorbed if the horizon exists. The reflected waves will be 
partly reflected back again at the angular momentum barrier but partly transmitted. This 
process occurs repeatedly, which phenomenon is dubbed as GW “echoes” [499, 500] (see 
also Ref. [501] for a simple explanation with the causal structure of the Green’s function 
of the wave equation with two potentials having disjoint bounded support). Therefore, the 
detection of echoes implies the existence of an exotic compact objects. The simplest model 
of the compact objects is that the horizon is replaced by a reflecting membrane at ~ Planck 
proper length from the horizon radius in Kerr spacetime. Abedi et al. have analyzed three 
O1 BBH events and reported that they found echo signals at 3c significance level (0.011 in 
p-value) [502]. In their study, they assume the GW waveform in the frequency domain as 
N 
AG) teh o ye ee eee (61) 
n=1 

where y is the overall reflection rate, Atecho is the interval between neighboring echoes, 
N is the number of echoes, and ¢(f) is the phase shift. The function ho(f) is a merger- 
ringdown waveform extracted by cutting off the inspiral part smoothly from the inspiral- 
merger-ringdown waveform. The parameters y and Atecho are the most relevant parameters 
to characterize the echo waveforms on which the significance of the signal depends most. 
Since Abedi et al. used only 32 seconds for the background estimation, Westerweck et al. [503] 
have improved the background estimation using 4096 seconds and shown lower significance, 
0.032 in p-value, with the same template used in the paper by Abedi et al. [502]. The 
Bayesian analysis is also performed with the same template and small Bayes factors are 
reported [504, 505]. The spacetime considered in Abedi et al. is exactly Kerr spacetime 
outside the surface, while the assumption of reflection rate in their study is not physically 
reasonable. Nakano et al. have derived the reflection rate at the angular momentum potential 
barrier numerically, and calculated a fitting formula for the reflection rate R(f) [506]. The 
reflection rate becomes frequency dependent, unlike in the assumption in Abedi et al., and 
behaves as a high pass filter, i.e., R(f) — 0 [R(f) — 1] for larger (smaller) f. As a result, 
the lower frequency part of the waveform is strongly suppressed in the template given by 
Nakano et al. compared to that given by Abedi et al.. The modified template becomes 


N 
ESTER Re) Y RUE) re er ee RUDI D) (62) 
n=1 


In this setup, Uchikata et al. have shown that the significance of echo signals becomes much 
lower, p-value is larger than 0.9 for eight BBH events in O1 and O2 [507], although neither 
reduction nor increase of the significance was observed by increasing the number of samples 
when the original waveform proposed by Abedi et al. is employed. 

So far, we have focused on a model dependent search. A morphology independent search 
was done by Tsang et al. [508, 509], where they assumed that echo waveforms for any 
model can be effectively described as superpositions of sine-Gaussian waveforms. In their 
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study, no significant echo signal is found in O1 and O2 events, including the BNS merger 
event. The reflection rate at the surface can also vary due to the model of the object. 
In the analytical template proposed by Refs. [510, 511], the reflection rate at the object’s 
surface is also assumed as a parameter. The reflection rate at the potential barrier is given 
from BH perturbations. Further details of testing exotic compact objects includes the echo 
phenomenon are reviewed in Ref. [512]. 


5.2. Dark matter 


Arguably, one of the most important questions in science is to elucidate the nature of dark 
matter. Identifying the nature of dark matter definitely brings huge impact on cosmology, 
particle physics, and astronomy. Many dark matter candidates have been proposed in the 
literature. Since different dark matter candidates yield different type of observational sig- 
nals, various experiments targeted for particular dark matter have been conducted. In this 
subsection, we give a brief overview of how the GW experiments can shed light on some 
dark matter candidates. 


5.2.1. Axion. An axion is a strong candidate for the dark matter[513]. To probe the axion 
dark matter, it is important to study electromagnetic waves [514] and GWs propagating in 
the axion dark matter [515, 516]. In particular, the gravitational Chern-Simons term allows 
us to probe the parity violation in the gravity sector. We consider the model specified by 
the Lagrangian (33) with V = m?¢?/2. 

It is instructive to write down the quadratic action deduced from Eq. (57) using the circular 
polarization modes as 

2 
so = 1e f Bk dt (: = ees) [thay — hin, . (63) 
4 Mp 

To avoid a ghost, due to the signature flip of the factor (1 — kaal $/Mp), we need the 
cutoff scale below which we can use this effective action 


8 2 3 
k < kg = 6.7 x 10°Hz (= =) Jd IE (64) 


The axion is coherently oscillating as $ = $9 cos(mt) with the amplitude determined by 


10719eV p 
c 2.1 x 107eV NI 65 
o POS ( m ) 0.3GeV /cm? (68) 


Here, we used the observed local dark matter density p. The coherent oscillation induces the 


parametric resonance of GWs with a frequency corresponding to the axion mass. Using the 
current upperbound £ < 108km [517], we can estimate the length scale R19 at which the 
GW grows ten times larger as follows 


10-19eV4? /108km\? /0.3GeV /cm? 
R10 — 5.2x 1078pc ( — ) ( ? =) - icu $ (66) 


Since the Jeans length rj of the axion dark matter, which is the coherence length of the 


axion oscillation 
= m z^ p Sa 
= 4.3 x 1073 (v) UU 67 
" ^ PE (I0 eV 0.3 GeV /cm? (67) 
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is much larger than the growth length, there occurs a huge enhancement of the amplitude 
of GWs, which gives rise to the strong constraint on the Chern-Simons coupling constant or 
the abundance of the axion dark matter. We refer the reader to the original paper [515] for 
the details. For the detection strategy to probe this effect is under investigation. 


5.2.2. Primordial black holes. Primordial black holes (PBHs) are a dark matter candidate 
that is not an elementary particle. The idea of PBHs dates back to 1971 when Hawking 
proposed a possibility that BHs are produced directly by the gravitational collapse of an 
overdense region in the primordial universe much prior to the recombination era. PBHs 
heavier than about 10!°g do not lose their mass significantly due to the Hawking radiation 
over the age of the Universe. Since such PBHs are non-relativistic objects, interact with 
baryonic matter only gravitationally, and do not emit light, they have been considered as a 
dark matter candidate. Since the original proposal of PBHs, various observational searches 
for PBHs have been continuously conducted for a wide PBH mass range. So far, there are no 
observational evidence of PBHs, and upper limits on the PBH abundance have been placed, 
which, for some mass range, excluded the possibility that PBHs comprise all dark matter 
under the assumption that all PBHs have the same mass. Yet, there is still an open window 
where PBHs can provide all dark matter (see Ref. [518] for the most recent summary). 

LIGO/Virgo detection of binary BHs provoked explosive research activities on PBHs (see 
Ref. [8] for a review in this respect). An exciting possibility is that these BHs (or a part of 
them) are PBHs. Indeed, features of the LIGO BHs such as comparatively large masses and 
small spins do not contradict with PBHs. It is known that a dominant formation channel 
of the PBH binaries is by a tidal force by the surrounding PBHs in the radiation domi- 
nated era [519]. This channel can explain the merger rate of the BH binaries derived by 
the LIGO/Virgo observations if the PBHs constitute about 107? fraction of dark matter. 
Alternatively, a conservative assumption that LIGO BHs were produced by the standard 
astrophysical processes yields an upper limit on the abundance of stellar mass PBHs. Irre- 
spective of whether LIGO BHs are PBHs or not, that PBHs can comprise at most ~ 1073 
of dark matter is the most severe constraint for PBHs in the stellar mass range. This vividly 
demonstrates the powerfulness of the GW observations to test PBHs as dark matter. It is 
worth to mention that LIGO/Virgo Collaboration extended searches of the PBH mergers to 
sub-solar mass range [520]. 

Although stellar mass PBHs are unlikely to explain all the dark matter, a question still 
remains that such PBHs constitute a tiny fraction of dark matter. Answering to that question 
definitely gives a profound consequence to cosmology. To this end, it is essential to discrim- 
inate between PBHs and the astrophysical BHs from GW observations. In what follows, we 
give a short review of the several ideas along this direction proposed in the literature. 


Mergers at high redshifts. Contrary to astrophysical BH binaries which formed long after 
the Big Bang, PBH binaries exist since the Universe was still dominated by radiation. Due 
to a broad distribution of orbital radius and eccentricity of the PBH binaries, merger rate of 
the PBH binaries continuously extends to high redshift even beyond z ~ 100. On the other 
hand, the merger rate of the astrophysical BHs sharply drops beyond z ~ 10. Thus, searches 
for mergers beyond z ~ 10 provide a clean test of the existence of PBHs. Reference [521] 
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has shown that the proposed space-interferometer called B-DECIGO is able to achieve this 
goal. 


Mass distribution. Mass distribution is another observable that can be used to test PBHs. 
Let mı and mz be masses of the individual BHs in a binary and R(m1, M2,t) be the dis- 
tribution of the merger rate of the BH binaries at cosmic time t. Notice that R includes all 
the mergers and the observational bias, which depends on the detectors, should be taken 
into account in translating R to directly measurable merger distribution. Since this bias is 
possible to handle, we here regard R as an observable quantity. 

Reference [522] showed that PBH binaries formed in the radiation dominated era gives a 
characteristic mass dependence to R as 


R = Cy(m1i)Y(m2)(m1 + ma)*. (68) 


Here C is a constant independent of m1, m», (m) is a function which depends on the PBH 
mass function, and o ~ 1. Quite interestingly, R almost linearly depends on the total mass 
mı + m» irrespective of the PBH mass function which differs for each inflation model and 
is unknown?. It is worth mentioning that other physical mechanisms leading to BH mergers 
predict different a values. For instance, the PBH binaries formed by two-body encounters 
in the low-redshift universe yield a ~ 1.43 [524]. Astrophysical BH binaries formed in dense 
star clusters give o ~ 4 [525]. BH binaries formed in mass-segregated environments such as 
galactic nuclei give a that vary with the total binary mass ms [526]. Thus measurement of 
a provides a unique test of the PBH scenario which is independent of the assumption for 
the PBH mass function. 

Alternative approach is to directly compare the observed mass distribution with the theo- 
retical ones in the PBH hypothesis computed for several representative PBH mass functions. 
'This program has been already applied to the observational data, and some PBH functions 
are excluded [527]. 


Spin of PBHs. Another potentially powerful GW observable to test the PBH scenario is 
spin distribution. GW observations are primarily sensitive to a quantity defined by Eq. (2). 
PBHs form binaries dynamically and there is no correlation between the spin of individual 
PBHs and the orbital angular momentum. Thus, probability density of veg is an even func- 
tion, which is a robust prediction of the PBH scenario. Concrete shape of the probability 
density of eg is determined by the spin distribution of the individual PBHs W (J). Formally, 
the spin distribution of PBHs at cosmic time ¢ is given by 

W(J4) = i dJ! QU, J',t) f PET (69) 

5en(J’) 

where P(ôm, J) is the probability distribution of the density contrast ôm and the angular 
momentum J of an overdense region that collapses to PBH if ôm > ên. Q(J, J’, t) represents 
the evolution of the PBH spin from its initial value J’ to J at time t. P(m, J) is deter- 
mined once the statistical properties of the primordial density perturbations are specified. 


? Recent study generalized the original work and found that a takes value different from unity if 
the dark matter perturbation provides the dominant torque to produce the angular momentum of 
the PBH binaries [523]. 
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All these quantities must be determined in order to derive W(J,t). In Refs. [528-530], anal- 
yses related to P(ôm, J) and Q(J, J',t) have been performed. For the Gaussian primordial 
density perturbation, the angular momentum of PBHs coming from P(ôm, J) is estimated 
to be at most a few percent [529, 530]. In Ref. [531], din(J) for PBHs formed in the radiation 
dominated epoch was obtained as 


Oth © 0.62 + 0.015a7. , (70) 


where ax is the Kerr parameter. This result shows that the threshold increases as the angular 
momentum is increased and the formation of highly spinning PBHs are suppressed com- 
pared to slowly spinning PBHs. On the other hand, PBHs formed in the matter-dominated 
phase, which could happen prior to the radiation dominated epoch, likely to have significant 
spin [532]. 


So far, we have considered GWs from PBHs in the LIGO frequency band. This band 
corresponds to the merger phase of the stellar mass PBH binaries. However, GWs in differ- 
ent frequency bands could be also generated. The PBH binaries in the inspiral phase, due 
to their enormous number and weak signals, constitute the low-frequency stochastic GW 
background. Furthermore, primordial density perturbations as a seed of the PBHs provide 
another stochastic GW background by the non-linear mode couplings. These low-frequency 
GWs are not probed by the LIGO-like detectors, but could be detected by future detectors 
such as LISA. 

In what follows, we introduce possible constraints on the abundance of PBHs expected in 
future experiments by evaluating energy-density spectra of stochastic GW backgrounds in 
two independent ways. We show that the experiments are sensitive to constrain the fraction 
to cold dark matter (CDM) for 107? < fpgu € 1 (107? < fppu € 1) in case of the GWs 
from coalescing events (curvature perturbations) [533]. 

By using the merger rate of PBH binaries [534], we can calculate the energy-density spec- 
trum of stochastic gravitational wave background (SGWG), which is discussed in Ref. [535] 
to be 


Mout _ | 
m E 
-f Rppu(z) d GW 6, qs. (71) 
0 


OQaw(v) = De (1+ z)H(z) dw, 


where drew (vs) is the energy spectrum of GW produced by a BBH coalescence [536, 537], 


Vs = (1 + z)v, and Veus is the cutoff frequency. 
On the other hand, according to Appendix D of Ref. [533], we obtain the dimensionless 
energy-density spectrum of the induced SGWB to be 


S (v= LI ) 4 ds aa (zs) 7 (72) 


2m Gx,s (Toq) 


x O0(2-—k), 

where v = k/2x denotes the frequency of GW, and the dimensionless wavenumber k—k / ko 
is introduced for simplicity. We used notations of physical quantities written in Ref. [533]. 

Assuming a null detection of the GW signals from coalescing events in Eq. (71) or curvature 
perturbations in Eq. (72), we can constrain the abundance of PBHs, which is plotted in 
Fig. 38. The solid lines give the expected upper limits on the abundance by the non-detection 
of the GW signals from the coalescing events, and the shaded regions mean the parameter 
space excluded by the non-detection of the GW signals from curvature perturbations. 
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Fig. 38 Expected constraints on the PBH abundance versus the PBH mass from the null 
detection of the two kinds of SGWB. The existing observational constraint (red dashed) is 
plotted for comparison [538]. This figure is taken by Fig. 7 of Ref. [533]. 


6. Conclusion 


Gravitational wave physics and astronomy have finally become science that can be com- 
pared with the actual observational data after a long period of preparation. In this paper, 
we introduced important findings such as the discovery of massive binary black holes and 
binary neutron stars, which bring us a new implication about physical properties of nuclear 
matter, heavy element synthesis in the Universe, and the origin of short gamma-ray bursts. 
It is expected that this trend will be further accelerated as the international network of 
gravitational wave detectors on the ground expands and space gravitational wave anten- 
nas extends. Further improvements in sensitivity are planned for ground-based gravitational 
wave interferometers, including the construction of the third generation detectors[433, 434]. 
An increase in sensitivity by an order of magnitude increases the chance to detect an event 
having a very high signal-to-noise ratio, and at the same time increases the event rate by 
a factor of 1000 or more, which enables highly statistically accurate analyses of gravita- 


tional waves. It is certain that future gravitational wave observations will revolutionize the 
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understanding of the physical universe that we have accumulated so far, powered by the 
improvement of the accuracy of gravitational wave source position and multi-messenger 
follow-up observations for gravitational wave events. It is clear that we have not seen the 
great potential of this new probe of physics, gravitational waves, yet. 

As the observation technology advances, it becomes more important to analyze and inter- 
pret the accumulated data correctly. The real benefit from gravitational wave observations 
can be extracted only when we can construct an appropriate physical picture to explain the 
observed signal, by combining the data from other observational means to follow-up gravita- 
tional wave detections. To that end, we need to further promote multi-messenger astronomy, 
and advance the theoretical research for understanding gravitational wave sources. In addi- 
tion, in order to detect gravitational wave signals that have not been discovered yet to 
open up a new paradigm, it is necessary to further develop gravitational wave data analysis 
methods, making full use of the accumulated theoretical knowledge. 
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